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MANAGEMENT AND COORDINATION 


The South Carolina NASA/EPSCoR program. Grant Number NCC5-643, transferred from the 
University of South Carolina (USC) to the South Carolina Research Authority (SCRA), effective April 1, 
2002. The SC EPSCoR State Committee and the SC Space Grant Consortium continued to implement the 
project management plan developed at the beginning of the grant (under NCC5-174) to provide 
consistent, efficient, and effective management of the NASA/EPSCoR program. 

The SC EPSCoR State Committee maintained overall responsibility for grant management and 
the SCRA assumed fiscal responsibility of the grant. In November 2002, the State Office of the SC 
EPSCoR program moved into the same building occupied by the corporate offices of the SCRA. SCRA 
is a 501(c)(3) organization whose mission is focused on technology-based economic development. The 
SC Space Grant Consortium continued to aid in the development of outreach and evaluation initiatives for 
the NASA/EPSCoR research projects. 

SC EPSCoR State Committee 


During the final reporting period, the SC EPSCoR State Committee, consisting of representatives 
from academia, state government, and industry, attended biannual meetings to stay informed on federal 
solicitations and reports and discuss South Carolina's EPSCoR initiatives. John R. Raymond, MD, Vice 
President for Academic Affairs & Provost at the Medical University of South Carolina, serves as Chair of 
the State Committee. 

Through an Executive Committee, the State Committee provided direction and oversight for the 
NASA/EPSCoR grant through monthly conference calls to Larry Druffel, PhD, Project Director, and the 
State Office of the SC EPSCoR program. The Committee assigned full management responsibility for 
the NASA/EPSCoR program to both Dr. Druffel and Mitchell Colgan, PhD, SC Space Grant Consortium 
Director. Scott Little, PhD, SC EPSCoR State Manager, and Tara Scozzaro, MPA, SC Space Grant 
Consortium Program Manager, assisted in the administration of the NASA/EPSCoR program. 

SC NASA/EPSCoR Program Management 

Effective April I, 2002, Dr. Druffel, State Director of the SC EPSCoR program, began serving as 
the SC NASA/EPSCoR Project Director. Dr. Druffel is responsible for implementing the terms and 
conditions of the cooperative agreement and overseeing the grant's research activities, reporting to the 
State Committee and NASA. Drs. Little and Druffel met weekly to discuss NASA/EPSCoR activities and 
specific achievements. 

NASAEPSCoR program initiatives in South Carolina, in coordination with the SC Space Grant 
Consortium, were centrally managed from the State Office of the SC EPSCoR program, which operates 
through four full-time staff members, including Dr. Little, State Manager, a Fiscal Manager, a Program 
Manager, and an administrative assistant. One undergraduate student provides part-time assistance in 
administrative activities. The State Office staff provided centralized fiscal and program monitoring, 
implementation of outreach activities, and serve as a liaison to faculty and student participants. 

The State Office worked closely with members of the state’s Commission on Higher Education 
(CHE), representatives from the PhD granting institutions, and state congressional leaders to secure state 
funds to help meet the cost-share commitments required on active EPSCoR awards in the state. During 
the final reporting period, more than S350K in state matching funds were awarded to the research 
activities supported under this NASA/EPSCoR grant. 

During the final reporting period, the State Office continued to maintain NASA/EPSCoR 
program data; updated the SC EPSCoR website; arranged logistics for meetings, conferences/workshops/ 
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symposia, site reviews, and presentations to the SC General Assembly; prepared annual reports; 
established, maintained and monitored all federal and cost-share accounts; and managed review logistics 
for the participant support programs. These programs include seed funding for multi-institutional 
sustainable research clusters conducting research in support of NASA's strategic research areas, while 
simultaneously supporting SC EPSCoR goals. The SC EPSCoR program and the SC Space Grant 
Consortium announce the annual call for applications in early spring, with receipt of applications by late 
spring in the Space Grant Consortium office. The State Office of the SC EPSCoR program identifies and 
confirms at least three faculty with expertise in the technical area of the proposal to conduct a thorough 
review of the application by mail. The State Office sits a three-member panel to resolve any conflicts 
from the mail reviews, rank the proposals, and make funding recommendations. 

SC Space Grant Consortium 

The SC Space Grant Consortium and the SC EPSCoR State Committee worked closely to 
organize, plan, and manage South Carolina's NASA/EPSCoR program. Dr. Colgan and Ms. Scozzaro 
worked with Drs. Druffel and Little to oversee the progress of the program and assist in the evaluation 
and assessment of the research projects. The SC Space Grant Consortium continued to help integrate the 
NASA/EPSCoR program activities into its outreach and educational activities. 

The SC Space Grant Consortium, through seminars and faculty exchanges, also disseminated 
information about EPSCoR-related research to institutions and scientists not presently involved in the 
NASA/EPSCoR program. Members of the SC EPSCoR State Committee and State Office were invited to 
attend SC Space Grant Consortium biannual meetings and were included on the consortium's mailing list 
to be kept abreast of the latest news and opportunities from NASA. 

Project Milestones and Reporting 

SC NASA/EPSCoR grant recipients adhered to several reporting requirements. In addition to 
submitting financial and technical progress reports, the SC Space Grant Consortium and the 
NASA/EPSCoR Directors prepared semi-annual progress reports for the SC EPSCoR State Committee. 
These reports included details on research tasks, projected milestones, resources to be used, and 
anticipated results, and were reviewed by the State Committee during their biannual meetings. In 
addition, the Principal Investigators (Pis) of the research activities were required to present a formal 
progress report to the State Committee on an annual basis and to provide periodic “bullets of 
achievement.” These reports and achievement data served as a forum to review progress and plan for 
future activities. 

During the final reporting period, Dr. Colgan continued to maintain a NASA/EPSCoR program 
file that serves as a reference tool for project staff, SC EPSCoR and Space Grant representatives, and 
NASA. This file contains copies of the grant instrument, copies of the technical and cost proposals, 
financial reports, progress reports, correspondence, meeting summaries, technical notes, and other related 
project information. 

NASA/EPSCoR Pis prepared quarterly cost summary reports, which were used by Drs. Druffel 
and Little to track overall program status and to meet NASA disbursement reporting requirements. Dr. 
Little and his staff also identified expenditure trends, potential problem areas and their effects on the 
project timetable and funding schedule through these reports. 

Subject Inventions: None 

Inventory Report of Federally-Owned Property: None 
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Trident Research Center 
5300 International Blvd. 
North Charleston, SC 29418 
TEL (843) 760-2700 
http://www.scra.org 


June 30, 2004 


Diane Detroye, Technical Officer 
NASA Headquarters, Code FE 
Washington, DC 20546-0001 



RE: South Carolina NASA EPSCoR Grant Number NCC5-643 Matching Funds Report 
Dear Ms. Detroye: 

Consistent with the Special Conditions of NASA Cooperative Agreement “South 
Carolina EPSCoR Grant” Number NCC5-643, the cost share commitment for the period 
4/1/2002 - 3/31/2004 is $1,790,576.00. 

To date, $1,791,056.79 non-federal funds have been matched to SC NASA EPSCoR 
Grant Number NCC5-643. 

If you have any questions or need any additional information, please do not hesitate to 
contact K. Allyn Graham at 843-760-3235 or me at 843-760-3369. 



Chief Financial and 
Administrative Officer 
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Development and Enhancement of 
Research Capability for Aircraft 
Structures and Materials 


Summary Report for Year 5 Activity 
University of South Carolina Cluster on Aeronautics 


Dr. Michael A. Sutton, PI 
Department of Mechanical Engineering 
University of South Carolina 
Columbia, SC 29208 


Narrative 

This is the final report for the USC portion of the 1996-2002 (1 year extension) South 
Carolina NASA EPSCoR program, with consortium agreement number NASA NCC-5-174. 
Even though recent changes within the NASA Langley Research center organization have 
resulted in major changes in focus for both our program and NASA Langley’s areas of interest, 
without any doubt the USC portion of the EPSCoR program has been a complete success. In the 
last year of the program, the following areas have been the focus of our research activity work; 

1 . Damage Tolerance studies for friction stir welded joints 
(Michael Sutton, Anthony Reynolds) 

2. Applications of 3D computer vision to structures (Michael Sutton, Jeffrey Helm, 
Stephen McNeill) 

3. Mixed mode fracture experiments (Michael Sutton, John Ward) 

4. Mixed mode fracture simulations through development of a computational 
predictive tool for general 3D crack simulations in aerospace structures (Michael 
Sutton, Xiaomin Deng, Jianzheng Zuo) 

5. Development of 3D measurement tool for micro-scale structures 
(Michael Sutton, Dorian Garcia, Hubert Schreier) 

6. Development of a 3D measurement tool for nano-scale structures 
(Michael Sutton, Hubert Schreier, Dorian Garcia) 

7. Development of nano-scale pattern application methods (Michael Sutton, Wally 
Scrivens) 

8. Development of molecular dynamic simulation capabilities (Xiaomin Deng) 

9. Ultra high speed computer vision for deformation measurements 
during impact (Stephen R. McNeill, Michael Sutton) 

Specifically, through the support provided by NASA EPSCoR, our faculty have developed 
national and international recognition in the following areas 

Dr. Anthony Reynolds Friction Stir Welding for Aerospace Structures 

Dr. Michael Sutton Computer vision for 2D and 3D deformation 

measurements in aero-structures 
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In addition, through the support provided by NASA EPSCoR, our faculty are beginning to 
develop national and international recognition in the following areas 


Dr. Xiaomin Deng 
Dr. Yuh Chao 
Dr. Stephen McNeill 
Dr. Michael Sutton 


Molecular Dynamic Simulations 
Weld modeling 
High speed computer vision 
Micro-scale and nano-scale measurement 
Technology 


Below is a summary of the papers, faculty and students supported by NASA EPSCoR 


Faculty: 


Staff: 

Post-doctoral fellows: 
MS Students 


PhD Students: 


A.P. Reynolds 

M.A. Sutton 

S.R. McNeill 

Wei. Zhao 

Xiaomin Deng 

J.S. Lyons 

M.A. Sutton 

Daniel Wilhelm 

Jianzheng Zuo 

Bangcheng Yang 

Ruiyun Wu (to graduate 8/2004 

Robert Taylor (graduated 2002 

Ning Yuan (graduated 2003) 

Darin Lockwood (graduated 2001) 
Elmoiz Mahgoub (graduated 2003) 
Hubert Schreier (graduated, 2003) 
Nicolas Comille (graduates, 12/ 2004) 
Paul Peterson (graduated, 2001) 
Weiming Lan ( in progress) 

Yunhui Yan (graduated 2004) 

Peng Cheng (graduated 2001) 


Scientific Accomplishments: 


1. Three-dimensional Stress and Deformation Fields Around Flat and Slant Cracks Under 
Remote Mode I Loading Conditions (Mahgoub, Deng and Sutton) 

The phenomenon of slant fracture observed in stable tearing tests of many ductile materials is 
investigated, where an initially flat crack, loaded under remote Mode I conditions, tends to grow 
into a slant crack and stay in the slant configuration until final fracture. In an effort to identify 
potential reasons why cracks prefer to grow in a slant manner, three-dimensional (3D) finite 
element analyses of crack-front stress and deformation fields in Arcan-type specimens 
containing a flat or slant crack are performed under elastic-plastic and remote Mode-I loading 
conditions. In particular, the crack-tip opening displacement (COD) at a position behind the 
crack tip, the mean stress, the effective stress, and a constraint factor (defined as the ratio of the 
mean stress and effective stress) are studied and compared for the two types of cracks. 
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Analysis results reveal several stress/deformation field variations around flat and slant 
cracks under identical remote loading conditions. First, close to the crack front, the COD of a 
slant crack is greater than that of a flat crack. Second, at the specimen’s mid-plane, a flat crack 
leads to a higher constraint value ahead of the crack than a slant crack. Third, the effective stress 
ahead of a slant crack is greater than that ahead of a flat crack, especially close to the crack front. 
The above results seem to suggest that slant fracture may be preferred because a slant crack 
enhances the driving force in the form of a higher near-tip COD value and because a shearing 
type of failure is promoted in the case of a slant crack compared to a tensile type of failure in the 
case of a flat crack 

2. A Novel Method to Investigate FRP Composite-Concrete Bond Toughness (Wan, 
Petrou, Harries, Sutton and Yang 

A novel experimental method using modified double cantilever beam specimens and a 
customized test frame are introduced to evaluate bond characteristics and toughness of fiber 
reinforced polymer (FRP) composite overlays and a concrete substrate under mixed modes 
loading. A computer vision system is used to measure the crack location, near-tip deformations 
and crack opening displacement during the crack growth process. Digital image correlation is 
used to determine the crack opening displacement (COD) for flaws growing in the vicinity of the 
FRP-concrete interface. 

Results from this study indicate that during crack growth, (a) the Mode I component of 
COD is dominant for all angles of specimen loading, (b) the magnitude of the local Mode I 
component of COD is maximized when good bond quality is present and crack extension occurs 
within the mortar/concrete near the FRP-concrete interface and (c) good agreement exists 
between independent energy release rate estimates based upon both an approximate elastic 
double cantilever beam formulation and also use of the measured components of COD in a 
classical linear elastic expression. 

3. Basic studies of ductile failure processes and implications for fracture prediction (Zuo, 
Deng and Sutton) 

Fracture of ductile materials has frequently been observed to result from the nucleation, 
growth and coalescence of microscopic voids. Experimental and analytical studies have shown 
that both the stress constraint factor and the effective plastic strain play a significant role in the 
ductile failure process. Experimental results also suggest that these two parameters are not 
independent of each other at failure initiation. In this study, a methodology for characterizing the 
effect of stress constraint (which is defined to be the ratio of the mean stress and the effective 
stress) on ductile failure is proposed. This methodology is based on experimental evidence that 
shows the effective plastic strain at failure initiation has a one-to-one relationship with stress 
constraint. 

Numerical analyses based on plane strain and three-dimensional unit-cell models have 
been carried out to investigate failure initiation of the unit cell under different constraint 
conditions. Results from the numerical studies indicate (a) for each void volume fraction, there 
exists a local failure locus in terms of mesoscopic quantities, a m and a e , that adequately predict 
incipient local micro-void link-up, (b) the results are fully consistent with a failure criterion that 
maximizes mesoscopic effective stress for a constant level of constraint defined by the ratio An,= 
Om/Oe, (c) for high to moderate constraint A m , the link-up envelope values for a m and a e are 
consistent with limit load conditions where the critical principal stress oi c corresponds to the 
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maximum principal stress in the loading history and (d) for low constraint, the link-up envelope 
values for cr m and a e correspond to link-up conditions having high levels of plastic strain and a 
principal stress oi that is lower than the maximum value for this loading history. 

Thus, the results suggest that a two-parameter ductile fracture criterion is plausible, such 
as critical crack opening displacement (COD) and constraint, for predicting the process of stable 
tearing in materials undergoing ductile void growth during the fracture process. 

4. Deformations in wide, center-notched, thin panels: Part I: Three-dimensional shape and 
deformation measurements by computer vision (Helm, Sutton, McNeill) 

The response of wide, thin, center-notched, 2024-T3 al uminum panels undergoing far-field 
tensile load is investigated. Three panels with a notch length to panel width of 0.33 and widths of 305 
mm, 610mm and 1016 mm were subjected to far-field tensile loading. As part of the experimental 
program, two pairs of cameras were configured into separate stereovision systems and used to capture 
simultaneously both the global response of the sheet and the local response near a notch tip. Global 
areas, ranging in size from 250 mm x 250 mm to 550 mm x 550 mm, were imaged for each panel. A 
second stereo-vison system recorded images of a small area, 10 mm x 20 mm, ahead of one notch tip. 
Post-processing of the stereovision measurement data from global and local systems using three- 
dimensional digital image correlation was used to obtain the complete displacement field at each point 
in the region of interest. 

In general, results demonstrate that the combination of stereovision and three-dimensional 
digital image correlation is capable of accurately measuring true, three-dimensional structural 
deformations in regions undergoing both large out-of-plane displacements and large displacement 
gradients. Furthermore, 3D measurements on the panel specimen near the grip location are shown to 
provide an independent assessment of the true boundary conditions, with specimen slippage clearly 
noted in the 1016mm specimen. 

Results from the extensive notched, wide panel experimental program demonstrate that (a) 
each panel has an initial shape that deviates up to 3mm from planarity, with the greatest deviations 
occurring at the center of the notch, (b) the global load-displacement response is essentially linear for 
load levels that are well beyond the onset of large, out-of-plane displacements in the notch region, (c) 
increasing the size of the notched, thin panel specimen results in distinctly different surface 
deformations and deformed shapes, with three separate maxima/minima in the out-of-plane 
component of the largest panel and (d) local strain measurements near the termination hole are 
consistent with a plane stress condition. The region where tensile opening strains are above 2% extend 
several millimeters ahead of the hole, while compressive strains parallel to the notch direction are 
contained within a few millimeters of the hole. The in-plane shear strains are concentrated along 
circular lobes at +/-45 0 from the horizontal direction, a trend which is generally consistent with plane 
stress conditions. 

5. Deformations in wide, center-notched, thin panels: Part H: Finite element analysis and 
comparison to experimental measurements (Helm, Sutton, McNeill) 

For comparison to the experimental measurements, finite element analyses were 
performed for each panel. Results indicate that (a) for the 305-mm and the 610-mm panels, the 
finite element predictions and experimental measurements for two of the displacement 
components (U(x,y), W(x,y)) are in good agreement, with the finite element predictions for the 
grip displacement (V(x,y=L/2) being 5% lower than the measurements. For the 1016-mm panel, 
the finite element predictions did not match the experimental measurements; the full-field 
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measurements demonstrated that the panel was slipped by 0.5 mm during the experiment, 
resulting in significant changes in the panel’s deformation. 

In summary, the results indicate that the FE method is effective in modeling the overall 
behavior of wide, flawed panels up to the onset of flaw growth. The combination of FEA and 
three-dimensional measurements for accurate prediction of displacement boundary conditions is 
shown to be a viable hybrid approach for predicting flawed panel response. 

6. Advances in stereo light microscopy (Schreier, Garcia, Sutton) 

The increasing research focus on small-scale mechanical systems has generated a need for 
deformation and strain measurement systems for micro-scale applications. Optical measurement 
systems, such as digital image correlation, present an obvious choice due to their non— contacting 
nature. However, the transfer of measurement technology developed for macro— scale 
applications to the micro-scale presents unique challenges due to the differences in the required 
high— magnification optics. This article illustrates the problems involved in calibrating a stereo 
microscope using traditional techniques and presents a novel methodology for acquiring 
accurate, three-dimensional surface shape and deformation data on small-scale specimens. 

Experimental results demonstrate that stereo microscope systems can be accurately and 
reliably calibrated using \Latin{a priori} distortion estimation techniques in combination with 
traditional bundle-adjustment. The unique feature of the presented methodology is that it does 
not require a precision calibration target but relies solely on point correspondences obtained by 
image correlation. A variety of experiments is presented to illustrate the measurement 
performance of a stereo-microscope system. It is shown that the surface strains obtained from the 
full-field, three-dimensional measurements on tensile specimens undergoing large rigid body 
motions are within $\pm50$ microstrain of strain gage measurements for strains ranging from 0 
to 2000 microstrain. 

7. Effect of Higher Order Displacement Fields on DIC Displacement Component 
Estimates (Schreier, Sutton) 

Digital image correlation techniques are commonly used to measure specimen displacements by 
finding correspondences between an image of the specimen in an undeformed or reference 
configuration and a second image under load. To establish correspondences between the two 
images, numerical techniques are used to locate an initially square image subset in a reference 
image within an image taken under load. During this process, shape functions of varying order 
can be applied to the initially square subset. Zero order shape functions permit the subset to 
translate rigidly, while first order shape functions represent an affine transform of the subset that 
permits a combination of translation, rotation, shear and normal strains 

In this article, the systematic errors that arise from the use of undermatched shape 
function, i.e., shape functions of lower order than the actual displacement field, are analyzed. It 
is shown that under certain conditions, the shape functions used can be approximated by a 
Savitzky-Golay low-pass filter applied to the displacement functions, permitting a convenient 
error analysis. Furthermore, this analysis is not limited to the displacements, but naturally 
extends to the higher-order terms included in the shape functions. This permits a direct analysis 
of the systematic strain errors associated with an undermatched shape function. Detailed 
numerical studies are presented for the case of a second-order displacement field and first and 
second order shape functions. Finally, the relation of this work to previously published studies is 
discussed.cross-correlation coefficient 
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8. On systematic errors in digital image correlation (Schreier, Braasch, Sutton) 

Digital image correlation as a tool for surface deformation measurements has found widespread 
use and acceptance in the field of experimental mechanics over the recent years. The method is 
known to reconstruct displacements with a sub-pixel accuracy depending on various factors such 
as image quality, noise, and the correlation algorithm chosen. However, the systematic errors of 
the method have so far been neglected. 

In this work, we address the systematic errors of the iterative spatial domain cross correlation 
algorithm caused by grey value interpolation. We investigate the position dependent bias in a 
numerical study and show that it can lead to apparent strains on the order of fifteen percent of the 
actual strain level. Furthermore, we present ways to reduce this bias to acceptable levels. 

9. Banded microstructure in 2024-T351 and 2524-T351 aluminum friction stir welds. Part 
I: Metallurgical studies (Yang, Yan, Sutton, Reynolds) 

Results from careful investigation of the banded microstructure observed on horizontal 
transverse cross-sections in all AA2024-T351 and AA2524 aluminum FSW joints indicate the 
presence of periodic variations in (a) the size of equiaxed grains, (b) micro-hardness and (c) 
concentration of base metal impurity particles (e.g., constituent particles) that correlate with the 
observed band spacing. The latter trend is more distinct in AA2024-T351, which has a higher 
volume fraction of constituent particles resulting in easily recognized particle rich regions on 
horizontal cross-sections near the mid-thickness of the joint and well-defined variations in 
hardness. In AA2524, the trends are more muted but clearly visible. 

Results from recent numerical simulations of the FSW process enable interpretation of 
the trends in grain size along the weld centerline in terms of the time-temperature cycle 
experienced by the material. Specifically, the AA2524 FSW joints having low power and high 
input energy (i.e., the slow FSW) exhibit micron-size grain structure across both bands. 
Conversely, the fast and medium FSW in AA2524 have higher maximum temperatures and a 
corresponding six-fold increase in grain size. 

10. Banded microstructure in 2024-T351 and 2524-T351 aluminum friction stir welds. Part 
II: Metallurgical studies (Sutton, Yang, Reynolds, Yan) 

A series of micro-mechanical experiments have been performed to quantify how the 
friction stir weld (FSW) process affects the material response within the periodic bands that have 
been shown to be a common feature of FSW joints. Micro-mechanical studies employed 
sectioning of small samples and micro-tensile testing using digital image correlation to quantify 
the local stress-strain variations in the banded region. 

Results indicate that the two types of bands in 2024-T351 and 2524-T351 aluminum 
FSW joints (a) have different hardening rates with the particle-rich bands having the higher 
strain hardening exponent (b) exhibit a periodic variation in micro-hardness across the bands and 
(c) do not change the initial yield stress of the material in the different bands. 
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11. Mixed Mode I/II Fracture of 2024-T3 Friction Stir Welds (Sutton, Reynolds, Yang, 
Taylor) 

For the first time, a series of mixed mode I/II fracture experiments have been performed 
on both base material and three families of friction stir welds (FSWs) in 7 mm thick, 2024-T3 
aluminum plate; the FSW joints are designated hot, medium and cold due to the level of nominal 
heat input during the joining process. 

Results from the fracture tests indicate that the measured critical COD at a fixed distance 
behind the crack tip properly correlates both load-crack extension response and microstructural 
fracture surface features for both the base metal and all friction stir welds, providing an accurate 
measure of toughness. In addition, the COD values indicate that transition from mode I to mode 
II dominant crack growth occurs at lower loading angles for hotter FSW joints, with a truly 
mixed mode I/II COD measured during crack growth in the medium FSW joint. Also, using 
results from recent detailed FSW metallographic studies, specific features in the fracture process 
are correlated to the resulting FSW microstructure. Finally, the observed ductile crack growth 
path in all three welds tends to leave the under-matched FSW weld region as the far-field applied 
shear loading is increased, with the medium FSW being the only case where the flaw remained 
within the FSW region for all combinations of shear and tensile loading 

12. Mode I Fracture and Microstructure for 2024-T3 Friction Stir Welds (Sutton, 
Reynolds, Yang, Taylor) 

Detailed microstructural studies and mode I fracture experiments have been performed on 
both base material and friction stir welds (FSWs) in 7 mm thick, 2024-T3 aluminum plate. 
Three sets of FSW process parameters, designated cold, medium and hot to describe the relative 
amount of heat input to the FSW region, have been investigated. 

Microstructural studies indicate that (a) the FSW nugget grain structure is relatively 
uniform in all welds, (b) grain size decreases by 25% when the global power input is increased 
125%, (c) a banded microstructure exists on horizontal cross-sections traversing the weld region, 
where the band wavelength corresponds to the tool advance per revolution, (d) particles in the 
banded microstructure have the same elemental content as the impurities in the base metal, 
implying that the weld process is responsible for the observed particle redistribution and 
microstructural banding, (e) all welds exhibit a decrease in particle size and reduced particle 
volume fraction on the advancing side of the weld nugget and an attendant increase in particle 
volume fraction on the retreating side of the weld nugget, (f) hardness minima are present in the 
HAZ region outside of the weld nugget on both the advancing side and retreating side for all 
welds and (g) the hot weld has the lowest overall weld hardness. 

Results from mode I fracture tests indicate that (a) macro-scale parameters (e.g., COD, 
yield stress) appear to be adequate for characterizing the observed fracture behavior in the under- 
matched FSW, (b) critical COD measured at a fixed distance behind the crack tip is a viable 
fracture parameter for FSW joints, with COD demonstrated as being capable of properly 
ordering the load-crack extension response for both the base metal and the welds, (c) die FSW 
joint has a through-thickness variation in toughness, (d) the observed ductile crack growth path 
remains in the softer nugget region of the FSW throughout the fracture process, though the 
microstructure variations appear to have an effect on the local fracture process and (e) the 
fracture path can be correlated with the locations of the hardness minima and, for the medium 
weld, with observed changes in the weld nugget due to FSW process 
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13. Full-field Speckle Pattern Image Correlation with B-Spline Deformation Function 
(Cheng, Sutton, McNeill and Schreier) 

A full-field speckle pattern image correlation method is presented that will determine 
directly the complete, two-dimensional deformation field during the image correlation process on 
digital images obtained using computer vision systems. 

In this work, a B-Spline function is used to represent the object deformation field 
throughout the entire image area. This is an improvement over subset-based image correlation 
methods by implicitly maintaining position and derivative continuity constraints among subsets 
up to a specified order. The control point variables within the B-Spline deformation function are 
optimized iteratively with Levenberg-Marquardt method to achieve minimum disparity between 
the predicted and actual deformed images. Results have shown that the proposed method is 
computationally efficient, accurate and robust. The general framework of this method can be 
applied to n-dimensional image correlation systems that solve for multi-dimension vector fields. 

14. A Study of Residual Stresses and Microstructure in 2024-T3 Aluminum Friction Stir 
Butt Welds (Sutton, Reynolds, Wang and Hubbard) 

Three-dimensional residual stress mapping in an aluminum 2024-T3 arcan specimen, butt- 
welded by friction stir technique at a low welding speed (), was performed by neutron 
diffraction. Results indicate that the residual stress distribution profiles across the weld region 
are asymmetric with respect to the weld centerline, with the largest gradients in residual normal 
stress occurring on the advancing side of the weld. Relative to Longitudinal stress is the largest 
tensile component with a maximum tensile stress (~100 MPa) at 7.5 mm from the weld 
centerline, and larger residual stresses are located close to the top surface. 

Transverse stress is almost a tensile component but has a lower level of stress (max. ~60 
MPa), except that four compressive residual stress regions are obtained at ±12.5 mm from the 
weld centerline near both the top and bottom surfaces. The normal stress component is in the 
range from -40 to 25 MPa, with compressive stress regions across the weld centerline within 1 5 
mm. A relieve of the residual stresses in the section closer to the edge of the arcan specimen was 
observed. 

15. Basic Studies of Welds in a Tank Car Steel: Residual Stress Measurements and Weld 
Characterization for TC-128B Steel Plate (Abdelmajid, Sutton, Wang and Hubbard) 

Welding is the primary joining process used in the construction and repair of railroad 
tank cars. Since both fatigue crack growth and fracture of tank car structures are often initiated in 
the vicinity of welds, the effects of residual stress should be considered. Thus, for the first time, 
the enclosed work reports neutron diffraction results for all six components of the three- 
dimensional residual stress field on a transverse weld cross section in a widely used railroad tank 
car steel, TC128-B. 

Results for a non-heat treated specimen indicate that (a) the longitudinal residual stress 
approaches 70% of the uniaxial yield stress, (b) the residual effective stress ranges from 250 Mpa 
to 450 Mpa in the weld region and (c) the residual shear stresses are of the same order as the 
smaller residual principal stresses in the weld region. Finally, initial fatigue crack growth 
predictions for weld flaws indicates that the residual stresses can reduce the fatigue life by a 
factor of two. 
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1. E. Mahgoub, X. Deng and M.A. Sutton, “Three-dimensional stress and deformation fields around 
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Mechanics, 70 2527-2542 (2003). 

2. B. Wan, M.F. Petrou, K.A. Harries, M.A. Sutton, and B. Yang, “A Novel Method to Investigate 
FRP Composite-Concrete Bond Toughness”, ASCE Journal for Engineering Mechanics (in 
press). 

3. JZ Zuo, MA. Sutton and X Deng, "Basic studies of ductile failure processes and implications for 
fracture prediction". Fatigue and Fracture of Engineering Materials (in press). 

4. Helm JD, Sutton MA and McNeill SR, "Deformations in wide, center-notched, thin panels: Part 
I: Three dimensional shape and deformation measurements by computer vision", Optical 
Engineering, 42 (5), 1293-1305 (2003) 

5. Helm JD, Sutton MA and McNeill SR, "Deformations in wide, center-notched, thin panels: Part 
II: Finite element analysis and comparison to experimental measurements", Optical Engineering, 
42 (5), 1306-1320 (2003) 

6. H. W. Schreier, D. Garcia and M.A. Sutton, “Advances in Stereo Light Microscopy”, 
Experimental Mechanics: An International Journal, 44 (3) 278-289 (2004) 

7. H.W. Schreier, M.A. Sutton, "Effect of Higher Order Displacement Fields on Digital Image 
Correlation Displacement Component Estimates", Experimental Mechanics: An International 
Journal, 42 (3) 303-311 (2002). 

8. BC Yang, JH Yan, MA Sutton and AP Reynolds, "Microstructure and banding in 2024-T351 and 
2524-T351 aluminum friction stir welding joints. Part I: metallurgical studies, accepted for 
publication in Material Science and Engineering A, 364, 55-65 (2004). 

9. MA Sutton, BC Yang, JH Yan and AP Reynolds, "Microstructure and banding in 2024-T35 1 and 
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10. M.A. Sutton, A.P. Reynolds, B. Yang and R. Taylor, "Mixed mode I/n crack growth in 2024-T3 
aluminum friction stir welds". Engineering Fracture Mechanics, 70, 2215-2234 (2003). 

1 1. M.A. Sutton, A.P. Reynolds, B. Yang and R.K. Taylor HI, "Mode I Fracture and Microstructure 
for 2024-T3 Friction Stir Welds", Material Science and Engineering A, 354, 6-16 (2003). 

12. M.A. Sutton, A.P. Reynolds, B. Wang and R. Taylor, "Microstructure evolution in 2024-T3 
aluminum friction stir welds and effect on mixed mode fracture", Material Science and 
Engineering A, Vol A323, 160-166 (2002). 

13. M.A. Sutton, A.P. Reynolds, D.Q. Wang and C. Hubbard, "A Study of Residual Stresses- and 
Microstructure in 2024-T3 Aluminum Friction Stir Butt Welds", ASME Journal of Engineering 
Materials and Technology, Vol 124, 215-221 (2002). 

14. P. Cheng, M.A. Sutton, S.R. McNeill and H.W. Schreier, ’’Full-field Speckle Pattern Image 
Correlation with B-Spline Deformation Function", International Journal of Experimental 
Mechanics, 42(3) 344-353 (2002). 

15. M.A. Sutton, W. Zhao, C. Hubbard and D.Q. Wang, " Basic Studies of Welds in a Tank Car 
Steel: Residual Stress Measurements and Weld Characterization for TC-128B Steel Plate, ASME 
Journal for Pressure Vessels and Piping, 124 405-414 (2002). 
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South Carolina NASA EPSCoR Program at Clemson University 
Report on Year 5 Activity 


The research being conducted at Clemson can be broadly characterized as development of advanced 
materials, analysis methods, and design concepts for improving the performance of aircraft structures. 
Material investigations include studies of advanced titanium alloys, TiC reinforced titanium composites, and 
hybrid mixed-tow epoxy composites. Analysis methods are being developed that are applicable to these and 
conventional advanced materials. The effect of cyclic loads and presence of initial damage and defects on 
material behavior and structural response, as well as the damage tolerance of these materials and structures 
subject to monotonic and cyclic loads are common elements in many of the investigations of the four faculty 
participants and their graduate students. The effort is approximately equally divided between experimental 
and analytical investigations. A description of progress made during the past year on the four research 
projects is given below. 

(1) Damage Progression Analysis, Residual Strength Prediction, and Finite 
Element Development for Laminated Plates 

Sherrill Biggers, Jr., Professor of Mechanical Engineering, Clemson Univ. 

De Xie, Ph.D., Engineering Mechanics, currently Univ. Michigan 
Mohan Rathinasabapathy, M.S., Mechanical Engineering, currently ABAQUS Inc. 


Summary 

The objective of this project is to develop efficient and accurate analysis methods that can be used to 
predict the progression of damage in laminated composite structures and the residual strength of damaged 
structures. In particular, the interaction of local damage and global structural response is of primary interest. 
Application of these analysis methods to develop design concepts that improve damage tolerance and 
residual strength is also a goal of this project. Therefore, both computational efficiency and accuracy are of 
major concern. Damage can progress in different ways in different plies according to the ply properties and 
fiber orientation, the ply-level stress state and the condition of the interlaminar bond. Questions arise as to 
how damage affects the local properties of the material, how these changes in properties affect the global 
response of the structure, and how the global response affects further development of damage. Models are 
needed which are both computationally inexpensive as well as accurate. Research in the first four years 
included development of finite element analysis methods for postbuckled composite plates using a higher- 
order zig-zag theory and investigation of damage progression in tailored plates with and without cutouts. 
During the past year, work has focused on (1) the efficient and accurate modeling of delamination - 
progression in plates and sandwich panels and (2) modeling of through-the-thickness stitching in a way that 
can be incorporated in the delamination modeling. Dr. Xie was supported by this grant and received his 
Ph.D. in December 2002. Mr. Rathinasabapathy was also supported by this grant and received his Masters in 
August 2003. 

Section I: Narrative 

Description of Research and Accomplishments 

The research on delamination progression builds on investigations during previous years which 
focused on (1) evaluation of local damage evolution models, (2) simulation of damage progression under 
both monotonic and cyclic loading, (3) damage progression in postbuckled plates, and (4) damage 
progression and residual strength of tailored flat and curved plates with a cutout. During the past year 
particular attention has been given to (1) modeling of delamination growth after buckling, and (2) modeling 
stitching as a delamination growth control device. These last two topics are discussed briefly. 
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(1) STRAIN ENERGY RELEASE RATE CALCULATION FOR A MOVING 
DELAMINATION FRONT OF ARBITRARY SHAPE 
BASED ON THE VIRTUAL CRACK CLOSURE TECHNIQUE 


Introduction 

The onset and growth of delaminations are of great interest in the applications of composite 
structures since delaminations are one of the major failure modes. A delamination has usually been treated 
as a 2D or 3D crack so that facture mechanics can be utilized. Among the variables of fracture mechanics, 
the strain energy release rate (G/, G u and G H/ ) has been used most frequently as the criterion in determining 
the initiation and growth of delaminations in composite materials. Several closed-form expressions for (G) 
have been derived through analytical approaches for specific cases IM1 . For more general and complicated 
problems in engineering applications, numerical approaches such as finite element methods are required in 
conjunction with the virtual crack closure technique that was developed to calculate the strain energy release 
rate numerically. 

The virtual crack closure technique (VCCT) was proposed for 2D crack configurations by Rybicki 
and Kanninen [s * and was extended to three dimensions (3D-VCCT) by Shivakumar, Tan and Newman l6] . 
Usually VCCT is applied by using orthogonal meshes that have sides either normal or parallel to the 
delamination front. This can be easily achieved when the delamination front is fixed by adjusting FEA 
meshes to fit the pre-existing shape of the crack front, though this kind of matching between meshes and the 
geometry may require fine mesh patterns. 

When studying delamination growth in composite structures, using orthogonal mesh patterns 
becomes problematic due to the changing position and shape of the delamination front. This kind of mesh 
adjustment together with the associated nonlinear iterations makes the analysis prohibitively large in most 
cases. Furthermore, delaminations in practical situations are not likely to grow in a self-similar manner due 
to the uneven distribution of the strain energy release rate (G) across the delamination front and due to the 
interaction of local and global response. For example, a straight delamination front does not stay totally 
straight in a double cantilever beam; a circular delamination tends to become an elliptical one under 
compression loading. A delamination due to impact damage is often irregular and grows in an irregular 
manner. It is not practical to create a new mesh for the newly predicted delamination front at each increment 
in order to keep the meshes orthogonal to the delamination front. This is a major obstacle in studying 
delamination growth in realistic cases. 

This study presents a method to calculate the strain energy release rate at the delamination front 
accurately and efficiently without using orthogonal meshes. Towards this end, an 18-node special element 
was developed to calculate the strain energy release rate based on 3D-VCCT once the delamination front lies 
within the element. 

Methodology 

3D-VCCT and G-calculation Interface Element 

The 3D-VCCT can be used to calculate the components of strain energy release rate (G/, Gu and G m ) 
for node N shown in Figure 1 with an orthogonal mesh are as follows: 

G " (1) 
g, “ 

where cr 2 (r, s ) , (r, s ) and t s . (r, s ) are the stress components ahead of the crack front and 

w(A 1 ,s),u r (A 1 ,.s) and t/^ApS) are the components of crack openings behind the crack front. The 
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parameters r and s are the local coordinates normal to and along the crack front, respectively. AA is the 
virtually closed area. 



Figure 1. G Calculation by 3D-VCCT at Node N Fig. 2 . Element definition and node numbering 


To evaluate equations (1) numerically, a special 18-node element is proposed herein. Figure 2 shows 
the element definition and the node numbering. The element consists of two sets of node groups: the top set 
(filled circles in Figure 2) and the bottom set (open circles in Figure 2). When the element is placed at an 
interface between layers, the nodes in the top set are from the upper portion of the interface while the nodes 
in the bottom set are from the lower portion of the interface. The elements above or below an interface 
element can consist of four 8-node solids or four 4-node shell elements. If the interface is perfectly bonded, 
i.e., there is no delamination, the nodes in the top set will have the same coordinates as their corresponding 
partners in the bottom set. In this case, the gap between the top node set and the bottom node, set, 
exaggerated in Figure 2, does not exist. The element has zero thickness. 

Since the element has been defined to calculate the strain energy release rates Gi , G n , and G m at the 
point where node pair (1,2) is located, three very stiff springs are placed between nodes 1 and 2 to evaluate 
the forces in the springs in the x, y, and z directions in the global coordinate system. In all the applications of 
this element made to date, the value chosen for the spring stiffnesses has not been a critical choice. Large 
values should be selected and no numerical solution problems have been detected for a wide range of values. 
Thus the element is robust in its numerical application. Let k x , k y , and k : be the stiffness of the springs. The 
total strain energy in the spring is 

U = l k M -u 2 f + £*,(v, -v 2 ) 2 +\k z (w x -w 2 ) 2 

The first variation of U is 

SU = k x {u x -u 2 )(Su x -^ 2 ) + ^(v, -VjXcSv, -<Sv 2 ) + £ z (w, - w 2 )(Sw x - dw 2 ) 

= {&i}[K]{u} r 

where {<5u} = {5u x ,dv x ,5w x ,du 2 ,dv 2 , S vv 2 } , {u} = {u , , v, , w, , u 2 , v 2 , w 2 } . The stiffness matrix for the 
element is 



where 

k x 0 O' 

[K,]= 0 k y 0 
0 0 k z 

The stiffness matrix [K] is implemented within ABAQUS by means of a user defined element subroutine, 
UEL. 
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All the information necessary to calculate the three energy release rate quantities is present in the 
element without additional post processing. The further advantage of this element is its capability of 
predicting delamination growth since an energy release rate based criteria can also be easily implemented 
into the user subroutine of the special element. 


Calculation of G for a Delamination Front having an Arbitrary Shape 

The delamination front within an element that is centered on node N can be located by defining the 
two of the eight vectors pointing away from point N (Figure 3(a)) that divide the bonded nodes from the 
debonded nodes. Let these two vectors be defined as R A and R e where the “b” and “e” subscripts stand for 
the "beginning” and "end” of the bonded region where one moves counterclockwise from R 4 to R e in 
sweeping out the bonded region. The damage state for the surrounding nodes is defined by a variable 
referred to herein as the damage index, 

f 0 for perfectly bonded nodes 

m i ~ 1 

[1 for separated nodes 

The center node N is considered to be bonded. Knowing the damage indices, the delamination front can be 
located within an element having one or more values of damage index equal to one by defining the two 
special vectors 

8 8 

= 5><-,0 - R e = £(l - "iM 

i= 1 i=l 

Here, m 0 = m s and m 9 = m\. These two vectors form a close approximation to the delamination front that 
passes through point N within the element. Figure 3(b) gives one example in which all nodes except 3 and 4 
are bonded and therefore b = 5 and e =2. 



Figure 3. Vectors that Form a Delamination Front 

Determination of the vector that is normal to the delamination front at N is fundamental for the 
further calculation of strain energy release rates. The normal vector here is approximated as the resultant of 

the two vectors R A and R e . Following the right-hand-rule for the cross product, we have 

_» 8 8 

h b = kxR 6 = k x 5]w,_,(l - m,)(Rj + R iy }) = £/w,_,(l - m,)(-/y + /?J) 

i=l i=l 

8 8 

n e = -k x R e = -k x £(1 - + Rj) = £(1 - m f )m M ( - RJ) 

i=i i=i 

where the unit vectors i and j correspond to the natural coordinate (£, 77 ) directions and k is the unit 
vector normal to i and j . Therefore, the sum of the two vectors is 
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8 

J=1 

The vector is then normalized as a umt vector. The unit normal vector and its corresponding tangent vector 
are 


n i+n e 


P*+“. 


t = kxn = -n i + n j 

y u 

The three umt vectors (n,t,k) form a base for the local coordinate system. The components of strain energy 
release rate (Gi, Gn and G/jj) will be calculated within this local coordinate system. 

To determine Gj, Gu and Gm, the virtually closed area needs to be calculated first. The calculation is 
carried out in the space of the natural coordinates (£, 77 ) , see Figure 4(a). The shape functions for the 
element in Figure 4(b) are 

^2 (£> 7) = y (1 - £ 2 X- 1 + 77)>7; W 3 (£, 7 ) = t (1 - 77 2 )(- 1 + £)£ 

A7 4 (<f,77) = i(l-£ 2 )( 1 + 77)77; iV 5 (£,,7) = X( l-^ 2 )( l + £)£ 

N 6 {§ , 7 ) = |(1 - £)(1 - 7 )^ 7 ; N, (§, 7 ) = -|(1 + £)(1 - 7 )^ 7 ; 

N* (£ 7) = 1 (1 + £)0 + 7)^7; (£ 7 ) --1(1- <?)(! + 7)^7 



(a) (b) 

Figure 4. Scheme of Virtually Closed Crack Area 


Referring to Figure 4(a), the virtually closed area can be divided into two parts: A l and A 2 . Four 
points denoted as P* ( k = 1,2; j = 1,2, 3, 4) enclose each of the two sub-areas. The areas are functions of the 
coordinates of these points, i.e., A t = A t (P* , P 2 , P 3 * , P 4 * ), it = 1, 2 . 

For the area A t , the three known points are 

P ' = (0.°) Pi = ( 6 — 1 » 7e— , ) Pi = (!(£-, + & H ), l(7 e -l + 7e+l)) 

and the following three points are known for defining area A^ . 

P 2 = (0,0) P 2 2 = (T(^_j +£ i+ 1 ),l(7 A _, + 7*+i) P 3 2 = ( 6 + 1 , 74 + 1 ) 
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To obtain the two unknown points P 2 and P 4 2 , two temporary points P 2 and P 4 2 are obtained first. 
The point P 2 is the intersection of (1) the line parallel to the vector n passing through point P,'and (2) the 
line parallel to the vector R A passing the point P 3 2 . Therefore, the two linear equations 

n Ai -M2 1 =o 

R Jl ~ R Jll= R eyZ'-X ~ R exVe-X 

can be solved to obtain the coordinates of P 2 ' as 

. a i 

t = f( R eyte-X - R a rie-x) 

where 


d X = ~»x R ex +n y R ey- 

The same procedure is then applied to obtain the coordinates of P 4 2 as 

£4 = ~T ( R bv^b+\ ~ RbxVb+l ) 

d 2 ' 

~ ~f~( R by^b+l ~ Rbxtfb+X ) 


where 


d , =~n r R bx + nj? 


y* V by ■ 


As an approximation, the actual P 2 and P 4 are taken as a single point that is the average of points 


P' and P 4 , that is, 




Now all the points P* = (^* , 7 * ), (j = 1,2,3,4; A: = 1,2) necessary to calculate the virtually closed area 
have been determined in terms of the parameters (£, 77 ) . Using the following relations 

x) = N i ^ k j ,r 1 k j )x i y)=Ntf k j ,T 1 k j )y i 

points (£*, 77 *) can be transformed into points in physical space (x k ,y k ) where (3c,, y,) are the nodal 
coordinates in physical space. Once the points (x k ,y k )are determined, Gauss Quadrature is used to 
calculate the areas Ai and A 2 with the total virtually closed area A given by 

A = {c, A, (P, 1 , P' , P 3 ' , Pi ) + c 2 A 2 (P , 2 , P 2 2 , P 2 , P 4 2 )} 

where c/ = 1.0 when point P 4 does not correspond to a node or 0.5 when the point does correspond to a 
physical node. The same can be said for c 2 but the governing point is P 2 . These two points, as well as point 
P, 1 , are always points on the delamination front. The c, constants as defined here prevent double counting of 
the virtually closed area as the special element is centered at each node on the delamination front. 

Next, the three nodal force components at node pair (1,2) are computed by 

K = k A u i ~ u i) F y = k y (v 2 - vj F z = k z (w 2 - w,) 

F, is required to calculate G/ and can be used directly as calculated above. F x and F must be 
transformed into the local system to calculate G„ and G m as 

F 2 = n x F r + n„F„ 


y y 


F z = n y F x - n x F y 
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The opening behind the delamination front is defined as the point anti-symmetric to points Pj and 
P 4 about ( 0 , 0 ) in (£, 77 ) space as 

p „ = (£>,%) = (-{(£ + 1 ?) -W 2 + m 2 )) 

Therefore, the displacement components for the delamination opening in the global system are 
Aw = N, (£ 0 , tj 0 ) Aw,. Av = N, (£ 0 , tj 0 ) Av,. Aw = N l (£ 0 , tj 0 ) Aw,. 

The Aw is required to calculate G/ and can be used directly as calculated above. The Aw and Av 
must be transformed into the local system to calculate Gu and Gm as 

Aw 2 = n x Aw + n y Av Aw 3 = n y Au - « x Av 

The strain energy release rate based on equation (1) can be approximated as the product of the nodal 
forces at node pair ( 1 , 2 ) and the nodal displacement opening at point P 0 (£ 0 ,77 0 ) in the region enclosed 
by A x and A, . Therefore, the equivalent of equation (1) is 

G, * ~T;Aw G„ * ^ F 2 Aw 2 G m « ~ ^Aw 3 


Numerical Results 

The validity of the present approach was first assessed by four examples, each with a fixed 
delamination front. Two of these delaminations represent mode-I cracks while the other two are 
characterized as mixed-mode cracks. They are selected to ensure the correct implementation of the 3D- 
VCCT into the user subroutine UEL of ABAQUS. 

Several delamination growth problems were then studied using the present approach. The criterion 
for determining the initiation and growth of delamination based on strain energy release rate is selected as 


( 




G, 


V u /c ) 


\P 


\~ncj 


’hi 


V ^ me J 


>1 


where Ej is the delamination growth parameter. G, c , G uc , and G U ic are the critical values corresponding to 
mode I, mode II and mode III cracks, respectively, and are assumed to be constant during delamination 
growth in this study. The exponents a, p, and y were taken as 1 as suggested by Kutlu and Chang [7 l Once 
the E d > 1 , the delamination grows and the stiffness matrix of the interface element is set to zero, i.e., 


[K] = [0] . No difficulties have been encountered due to sudden reduction of the stiffness to zero as opposed 
to a more gradual reduction used by others to aid in numerical convergence. 


Validation Examples for Fixed Delamination Front 

In order to evaluate the accuracy of the present approach, the method was first applied to four 
benchmark problems: an end delamination in a double cantilever beam (DCB) test specimen, an embedded 
elliptical crack in a solid body, a through-the-width delamination in a buckled plate, and an embedded - 
circular delamination in a buckled laminated plate. The first two are of typical mode-I dominated cracks 
while the later two are mixed-mode cracks. There are closed-form solutions or FEA predictions available 
for comparisons. 

Figure 5 shows the results for the DCB. The analytical solution given by Kanninen [8] used classical 
beam on elastic foundation theory to model the one of the two layers and therefore imposed the assumption 
of a straight crack front due to the 2D configuration. The three dimensional ABAQUS standard element 
"C3D8I" was used to model the specimen. There was only one element in each layer across the width so that 
the delamination front was forced to be straight, consistent with Kanninen’s assumption. Ten elements were 
used through the thickness of the delaminated layers and 90 elements were used along the length with the 
delamination occurring at the mid-surface. The mode-I strain energy release rate plotted in the figure was 
normalized with respect to the value G/ 0 as shown. Kanninen’s results were given in terms of stress intensity 
factors were reproduced and converted to G/ to prepare these plots. The specimen has a width (into the 
paper) of “B” and the remaining parameters are defined in the subplot in Figure 5. The solid line shows the 
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g-G/Go g=G/Gio 


analytical solution and the values predicted by 3D-VCCT are shown as black dots. The difference between 
the analytical solution and the 3D-VCCT approximation is within 1.58%. The larger errors occur for the 
shortest crack lengths where the effects of transverse shear deformations, not accounted for in Ref. [8] but 
included in the current approach, are the greatest. 

Figure 6 shows the distribution of the normalized Gj along the crack front of an embedded elliptical 
crack in a solid rectangular prism. The delamination shape and geometric parameters are shown in the 
subplot of the figure. The tensile stress a was applied in the direction normal to the paper at a distance far 
from the embedded delamination. The proportions used in this example for the characteristic dimensions are 
W=T, c=2a, and c/W=0.2. The "C3D8" element was also used for the solid body and an orthogonal mesh 


was used in the region of the delamination. Results are shown only for 0° < O < 90° due to the symmetry. 
The black dots show the distribution of Gj calculated from the 3D-VCCT and solid line represents the 
analytical expression given in Ref [9], The G, values plotted were normalized the by G m of circular 
delamination with radius “a”. Accuracy s imil ar to the previous example was observed, the difference 
between 3D-VCCT approximations and analytical solutions being within 1.33%. 

Figure 7 compares the strain energy release rates of a buckled, through-the-width delamination given 
by the analytical approach developed by Chai, Babcock and Knauss [1] with that given by the present 3D- 



Figure 5. Comparison of Analytical Solution and 
3D-VCCT for DCB 


Figure 6. Comparison of Analytical Solution and 
3D-VCCT for EEC in Solid Body 



(a) h/T=0. 1 (b) h/T=0.2 

Figure 7. Comparison of Analytical Solution and Present Approximation for a 1-D Delamination 


- 22 - 





VCCT method. The black dots are the VCCT results. The energy release rate shown in the figure is the total 
energy contributed by both mode I and mode H because separate values were not computed in analytical 
approach. The total G is plotted against the delamination length ( [I/L ) for two different delamination depths 
(h/T= 0. 1 and 0.2) and for three different compressive load levels (eg /e L ). The present approach based on 
3D-VCCT shows excellent accuracy compared to the analytical solutions. 

The final validation example is reported in Figure 8. Whitcomb 1 ' 01 also performed a finite element 
analysis based on 3D-VCCT. The sublaminate thickness (h) and the laminate thickness (H) are 0.4 mm and 
4 mm, respectively. The radius of the delamination (a) is 15 mm and the half width of the laminate plate (W) 
is 50 mm. The load level corresponds to e x = -0.005. A “homogeneous quasi-isotropic” material with 
(±45/0/90) s stacking sequence was selected in the study and the following mechanical properties for a typical 
graphite/epoxy were used: 

E„ = 134 GPa G 23 = 3.43 GPa v 23 = 0.49 

En = E 3 3 = 10.2 GPa Gn= Gn = 5.52 GPa v 12 = Vi 3 = 0.3 

The comparison in Figure 8 shows good agreement between 3D-VCCT employed by Whitcomb and the 
present 3D-VCCT method, especially in regions where the Gj and Gjj values are maximum, i.e. where growth 
would occur. Larger relative differences occur in areas where the Gj and Gn values are small. This is likely 
due to differences in treatment of the contact in this area. Slight relaxation of the gap closure constraint very 
close to the delamination front caused much better agreement with Ref. [10] in the regions of small G while 
not affecting the excellent correlation in other regions around the circumference of the delamination front. 



Figure 8. Comparison of FEA Solution in Ref [9] and Present 
Approximation for an Embedded Delamination 


Each of the four validation examples has shown excellent correlation between analytical solution or 
previous FEA analysis and the present approach. This confirms the accuracy of the 3D-VCCT for fixed 
delamination fronts and leads to consideration of the delamination growth studies discussed below. 

Delamination growth in a double cantilever beam 

Now consider the growth of a delamination in a double cantilever beam shown in Figure 9. The 
specimen is subjected to an end loading in the form of displacements that are constant across the specimen 
width. In itial studies were performed on other DCB models with mesh patterns aligned with the specimen 
geometry and oriented at 45 degrees to the specimen length. The current method produced strain energy 
release rates that were essentially independent of the mesh orientation. Here the method is applied to a 
specimen for which test results are available from Ref. [1 1]. The specimen is 150 mm long (L), 20 mm wide 
( B ) and composed of two 1.98-mm-thick ( h ) plies of unidirectional materials. The initial delamination length 
(a) is 55 mm. The material is_graphite/epoxy with the following properties: 
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E„ = 150.0 GPa 
E 22 - E 33 = 11.0 GPa 
G 12 — Gjj=6.0 GPa 
Gjj = 3.7 GPa 
V/2 ~ Vj 2 = 0.25 
V 2 J = 0.45 
G/c = 0.268 N/mm 
G/ic = 1 -45 N/mm 



Figure 9. Geometry of Double Cantilever 


The two models that were used are shown in Figure 10. ABAQUS standard “C3D8I” brick elements 
were used. Both models have 120 elements along the span and one element through the thickness of each of 
the two layers in the specimen. Two mesh sizes were used across the width: one (as in Ref. [1 1]) and ten 
elements, i.e., B01 and B10 mesh patterns, respectively. This allowed investigation of the importance of 
nonuniformity of the strain energy release rate across the specimen width and the associated curvature of the 
delamination front. 



Figure 10. Two Finite Clement Mesh Patterns 

Plots of the end reaction force versus end deflection are shown in Figure 1 1 for the two models and 
also for the test data from Ref. [11]. The distribution of mode I strain energy release rate across the 
specimen width is shown in Figure 12 for the load level that creates the first growth in the delamination 
front. The B10 model shows a significant variation G, across the specimen width. For this material, the G/c 
value is remarkably close to the uniform G/ predicted by the simpler B01 model. This indicates that while 
the simpler model would accurately predict initiation of delamination growth of a pre-existing straight-front 
delamination, it might not do so well as the delamination grows and the front deviates from its straight form. 
This fact is illustrated by the much better correlation of the B10 results with the test data than for the B 1 
results. The first open circle on the plot corresponds to the initial delamination growth and if the front is 
artificially forced to retain the straight delamination front, the simple model considerably underestimates the 
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Figure 13. Delamination Front Shapes During Growth for BIO Mesh Pattern 

ultimate load. The more refined model does an excellent job in predicting the test results. Since the current 
3D-VCCT approach does not require the elements to be orthogonal to the delamination front, it does an 
excellent job in dealing with the change in delamination shape as growth proceeds. Figure 13 shows the 
progression of the delamination where the slight deviation from straightness of the front can be seen. These 
10 patterns correspond to the 10 load-deflection points noted with circles in Figure 1 1. 


Delamination Growth under Compressive Loading 

To further study the sensitivity of the approach to the mesh pattern and the interface element 
orientation in the more difficult case of delamination buckling and post-buckling, three cases studies were 
examined. Their geometry and dimensions are shown in Figure 6 and Table 1. The delaminations are 
embedded in flat laminated plates and have the initial shapes shown in gray color in Figure 14. The three 
delaminated plates have the same 
stacking sequence. [0/90/90/D/0], with 
“D” referring to “Delamination”. 

Therefore, the delamination is located at 
the interface between the top 0° and 90° 
plies. The mechanical properties are 
given in Table 2. The compressive load 
is applied horizontally (1 -direction) on 
the edge of the plate as a uniform 
displacement. To initiate the out-of- 
plane deformation, an imperfection of 

2% the plate thickness is prescribed in the center of the plate. The imperfection is imposed upward for the 
sublaminate above the delamination and downward for the sublaminate below the delamination. This 
imperfection does not significantly affect the post-buckling behavior of the sublaminates. Geometrically 
non-linear analysis (NLGEOM) was carried out to include large deformations in the post-buckling stage. 

Each ply of the laminated plate is modeled 



Figure 14. Geometrical Configurations of Delaminations 


using ABAQUS brick elements (C3D8) across the 1- 
2 plane and one element through the thickness. To 
employ the VCCT, the interface elements are placed 

Table 1. Geometry Dimension 


Table 2. Mechanical Properties (Carbon/Epoxy) 



t/L 

a/L 

b/L 

Case (1) 

0.02 

0.4 

N/A 

Case (2) 

0.0333 

0.6 

0.2 

Case (3) 

0.0333 

0.8 

N/A 


Stiffness, Poisson’s Ratio 

Critical G (N/mm) 

Ei.(GPa) 

109.34 

Gjc 

0.306 

E 22 , E„ (GPa) 

8.82 

^-hic 

0.632 

G 12 , G ,3 (GPa) 

4.32 

^IIIC 

0.817 

G 23 (GPa) 

3.20 



V 12, V 13 

0.342 



V 23 

0.520 
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between the sublaminates to hold the nodes together on opposite sides in the interlaminar region. To support 
the general loading, three springs are placed at each node to act in mutually perpendicular directions 
corresponding to the three global coordinates. The stiffness of the springs in the three directions is chosen to 
be the same value as “k”. Due to the large out-of-plane deformation in die post-buckling stage, contact 
phenomena may be relevant in the delaminated regions. In order to prevent penetrations between the two 
contacting surfaces, gap elements (GAPUNI) were used to only allow separation of the node pairs in the out- 
of-plane direction. 


Throueh-the- Width Delamination 

The first example involves a through-the-width delamination that is investigated to gain some basic 
understanding of mesh sensitivity in a simple one-dimensional growth case. The geometrical configuration 
is shown in Figure 14(a). Due to the symmetry of the geometry and loading conditions, only one half of the 
structure is modeled. To perform the convergence study, four meshes with varying density were used and 
are shown in Figure 15. There is only one element across the width so that the delamination is forced to 
remain straight when it propagates along the 1-direction (local direction). For this specific mesh pattern, the 
enforced front approach and the front tracing approach using the adaptive normal vector direction predict 
identical results. Therefore, only one end load (N x ) vs. nominal strain (U/L) is plotted for each mesh density 
in Figure 16. Note that the four plots corresponding to four mesh refinements have the same initial slope 
prior to the buckling point. To help see the convergence clearly, the ultimate loads taken from the four plots 
in Figure 16 are plotted against the number of elements in Figure 17. It can be seen that the ultimate load 
decreases rapidly as the meshes are refined and tends to approach a constant after the number of elements 
reaches a certain degree of refinement. This indicates that the solution has converged. 



Figure 15. Four Refined Meshes for the 
Through-the-width Delamination. 



U/L 

Figure 16. Avg. Running Load N, vs. Nominal Strain 
U/L 


Another concern is the importance of choosing a proper stiffness value for the spring of the interface 
element. In order to study this, a series of calculations were performed with varying spring stiffness. The 
ultimate load (N x ) u is plotted against the spring stiffness (k) normalized with respect to the longitudinal 
stiffness (Ei i) in Figure 18. For each mesh density, there is a valid range for stiffness values that produce the 
same ultimate load. The calculations will encounter convergence difficulties if the stiffness falls out of the 
valid range. The finer the mesh, the broader the valid range will be. As expected, the stiffness used has 
almost no influence on the ultimate load for a wide range of numerical values. In other words, the approach 
is not sensitive to the interface element stiffness as long as a reasonable value is chosen. For this example, 
log(k/E n )=6.0 or k=109.34><10 6 GPa, which is valid for all of the mesh densities, would be one of many such 
reasonable values. 
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Number of Elements (N e ) 
Figure 17. Mesh Convergence Study 



Octagonal Embedded Delamination 

Next, consider an embedded delamination having an octagonal shape shown in Figure 14(b). Due to 
the symmetry, only a quarter of the plate has been modeled. Two mesh patterns and two interface element 
orientations were adopted. Investigations were made for combinations of the mesh patterns and the interface 
element orientations to study the approach's sensitivities to them. Mesh pattern A (MA) and mesh pattern B 
(MB) are shown in Figures 19(a) and (b), respectively. Each mesh pattern has four meshes of varying 
densities. 





A25 A100 A400 A1600 

(a) Mesh Pattern A (MA) 




B 50 B200 B800 B3200 

(b) Mesh Pattern B (MB) 

Figure 19. Mesh Patterns and Interface Element Orientations for Octagonal Embedded Delamination 
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Only orientation A (OA) for the interface element was used together with mesh pattern A. This is 
the baseline problem and only the adaptive normal vector approach is used. Four plots of average running 
load (Nx) vs. nominal strain (U/L) are shown in Figure 20 for four different mesh densities. Some important 
features can be observed from the figures. First, the greater the number of elements, the larger the range of 
loading over which delamination growth can be predicted and followed. For a coarse mesh such as A25, the 
initial load (Nx), at which delamination starts to grow, coincides with the ultimate load (Nx) u . As a 
consequence, there is no delamination growth observed before the ultimate load. For a refined mesh such as 
A1600, significant delamination growth occurs between the initiation and ultimate loads. Second, the 
method converges with respect to mesh refinement for both delamination growth initiation load (Nx)j and the 
ultimate load (N x ) u . These loads decrease rapidly as the meshes becomes finer but tend to approach constants 
after the number of elements reaches a certain degree of refinement. 



U/L U/L 

Figure 20. Avg. Running Load N x vs. Nominal Strain U/L (MA with OA), Front Tracing Approach, 

Octagonal Embedded Delamination 

To study the approach's sensitivity to the approximation of the delamination front normal direction, 
two orientations (OA and OB) of the interface element were applied to mesh pattern B. The plots of average 
running load (N x ) vs. nominal strain (U/L) curves obtained from the present approach are shown in Figure 
21. Again, the finer the mesh, the larger the load range over which the delamination growth can be traced. It 
is clear that, regardless of the mesh refinement, the approach predicts two nearly identical curves with the 
two different orientations of the interface elements (OA and OB). Since the delamination starts to grow at 
the root of the top horizontal edge where the interface elements have the same orientation for OA and OB, 
the same growth initiation loads are predicted for mesh B3200. The ultimate loads have a relative difference 
of 1.64%. This excellent agreement shows that the approach is not sensitive to the orientation of the 
interface element. The approach works well in detecting the delamination front. 

Comparing Figure 20 and Figure 21, the approach's sensitivity to the mesh pattern is clearly seen. 
There exists a significant difference between the ultimate end load vs. nominal strain curves of two coarse 
meshes (A25 and B50) because the problems have not converged yet at these mesh densities. As the mesh is 
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refined to A1600 and B3200, the curves get closer. The relative difference between the two delamination 
growth onset loads for B3200 and A1600 is 2.99%. For the ultimate loads compared to A1600, the values 
for B3200 with OA and B3200 with OB have relative differences of 2.85% and 1.19%, respectively. The 
predictions from the two mesh patterns eventually converge to approximately the same constant. The 
conclusion can be drawn that the approach is not sensitive to the mesh pattern once the mesh is fine enough. 



6.000 0.002 0.004 0.006 0.008 0.010 O.OOO 0.002 0.004 0.006 0.008 0.010 


U/L U/L 



6.000 0.002 0.004 0.006 0.008 0.010 6.000 0.002 0.004 0.006 0.008 0.010 


U/L U/L 

Figure 21. Avg. Running Load N, vs. Nominal Strain U/L (Mb with OA and OB), Front Tracing Approach, 

Octagonal Embedded Delamination 

It is desired to determine how the results would change if the normal vectors are forced to align with 
the mesh pattern established according to the delamination shape prior to loading. This approach is the 
“enforced front” approach. By comparing the dashed and solid lines in Figure 22, a significant difference 
can be seen between the results for the two different element orientation schemes. Furthermore, this 
difference is eliminated through mesh refinement. The enforced front algorithm obviously depends on the 
orientation of the interface elements. Therefore, assuming normal vectors prior to the analysis may prevent 
the correct result from being obtained even with very refined meshes. 

To summarize the above discussion, the ultimate loads (N x ) u are plotted versus the number of 
elements (N e ) in Figure 23. One can see that all three meshes (MA with OA, MB with OA, and MB with 
OB) predict very similar results through using the front tracing approach. If the front is enforced in an 
appropriate way (such as MB with OB), it may predict a good result. If the front is not enforced correctly 
(such as MB with OA), it may predict a poor result. Of course, if the front is enforced, one cannot be certain 
that orientations that are good ones prior to delamination growth remain as good ones after growth. To help 
understand this behavior, the delamination shapes at U/L=0.004 for different combinations of mesh patterns 
and interface element orientations are shown in Figure 24. Three combinations using the front tracing 
approach and MB with OB generate a similar shape for the delamination front at this instant, see Figures 
16(a)(b)(c) and (e). However, MB with OA using the enforced front approach predicts a totally different 
delamination shape as can be seen in Figure 16(d). Based on the above studies, the conclusion can be drawn 
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that with front tracing, all mesh patterns and interface element orientations, even the very simple mesh 
pattern A with “OA” oriented interface element can achieve very good results. 



U/L U/L 

Figure 22. Avg. Running Load N, vs. Nominal Strain U/L (MB with OA and OB), Enforced Front Approach, 

Octagonal Embedded Delamination 



Number of Elements (N e ) XI 0 3 Log(N e ) 

Figure 23. Summary for Octagonal Embedded Delamination 
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(a) A400, OA (b) B800, OA tel B800. OB 



(d) B800, OA (e) B800, OB 



Figure 23. Front shape for octagonal embedded delamination at U/L=0.004 

3.3 Square Embedded Delamination 

Finally, the same analysis procedure was also used to study a square shaped embedded delamination 
having the geometry shown in Figure 14(c). The definitions of the mesh pattern (MA and MB) and the 
interface element orientation (OA and OB) are the same as those for the octagonal embedded delamination 
discussed in the previous section. Again, "OA" oriented interface elements were imposed on mesh pattern A 
to provide a baseline problem and the front tracing approach was used. Results very similar to those 
discussed above were obtained for this delamination shape. As a summary, for all the cases considered here, 
the ultimate loads (N x ) u are plotted versus the number of elements (N e ) in Figure 24. These results once 
again demonstrate the ability of the current approach to capture significant aspects of the delamination front 
shape and location during delamination growth. 



Number of Elements (N e ) X10 3 Log(N e ) 

Figure 24. Summary for Square Embedded Delamination 
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Delamination growth in an irregularly shaped 2-D delamination 

As a more difficult and realistic test of the present 3D-VCCT method, a delamination with an 
initially complex shape and its growth to laminate failure are considered here. De Moura [!2! et al. presented a 
study of growth of a delamination created by a low-velocity impact event in a panel that is loaded in 
compression. Their experimental test data are used here for comparison. They also developed a 
computational method to track delamination growth based on an indirect usage of fracture mechanics, a 
material softening approach to damage development, and a penalty approach to prevent interpenetration of 
delaminated layers. The current approach differs in its direct application of mixed mode fracture mechanics 
through the 3-D VCCT method described earlier, immediate failure once the fracture criterion is satisfied, 
and standard gap elements to directly eliminate layer interpenetration. 



Figure 25. Delaminated Plate and Finite Element Mesh 

This application uses the same [O^OjJs graphite-epoxy laminate (see Figure 25(a)) loaded in 
compression in the 0°-direction that was presented in Ref. [12]. The laminate was impacted with a low- 
velocity mass resulting in the delamination between the outer 0/90 interface in the shape shown in Figure 
25(b). The linear combination of energy release rate ratios 

£j = C L . + c^ + 

G,c G IIC G mc 

was used as the fracture criterion. The same material properties, layer thicknesses, and initial imperfection 
stated in Reference [12] were used herein. The model used here contained 4 solid elements (ABAQUS 
C3DSI) through the thickness and a 20 x 20 mesh of elements across the plan of 'A of the plate. Symmetry 



End Displacement, U-, (mm) 

Figure 26. End Displacement vs. End Reaction 
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Figure 27. Contours of the U 3 During Delamination Growth 

conditions were used at the plate center lines and the unloaded edge is free. Interface elements were included 
at the upper 0/90 interface. 

Predicted load displacement data are plotted in Figure 26. The open circles are particular load levels 
selected for presentation of information on growth of the delamination in two later figures. The predicted 
ultimate load is within 1.3% of the measured test load reported in Reference [12], This is somewhat better 
than the correlation (3.8% error) with the method presented in that reference. This excellent correlation is 
remarkable considering the complex shape of the initial delamination and the complex pattern of growth. 

Deformed laminates are shown in Figure 27 for the 1 1 particular load levels identified in Figure 26 
as noted above. These plots are oriented in the same way as shown in Figure 14 with the load direction being 
from lower left to upper right. Plot (1) shows the delamination buckling after the delamination grows to 
include one additional separated node pair. It is interesting to note that the lower, thicker sublaminate also 
tends to buckle upward initially. The delamination grows first in the load direction and then generally in the 
direction transverse to the load direction as the delamination buckle increases in size and separation height. 
At some point close to ultimate load, the lower sublaminate buckles downward. Plot (9) corresponds to the 
deformed shape at ultimate load. When additional end displacement is applied, the delamination grows 
rapidly to the free, unloaded edge of the plate (plot (10)) and then continues to grow' until it encompasses the 
whole interface (plot (1 1)). 

The progression of the delamination can be seen in more detail in Figure 28. The 1 1 plots here 
correspond to the same 1 1 plots shown in Figure 27. Here the initially delaminated node pairs are shown as 
black circles and the shape of the initial delamination is shown. The open circles indicate delaminated node 
pairs at increasing load levels. The progression of the delamination is easy to follow in this set of plots. The 
progression predicted here is slightly different from the one plot reported in Reference [12] although the 
general trend and the global characteristics are the same. 


- 33 - 






Conclusions 

This section of the report has shown that the proposed algorithm for tracing the front of a 
delamination growing in a non-self-similar way is not sensitive to the values used for the interfacial spring 
stiffness, the orientation of the interface element, or the mesh pattern if the mesh is reasonably fine. To study 
these sensitivities, three cases having different initial delamination shapes were examined. Using a simple 
through-the width delamination, the spring stiffness value was shown to have almost no influence on the 
predictions as long as the numerical value is selected with a very wide range. For the octagonal embedded 
delamination, the ultimate loads obtained from different interface element orientations and the same mesh 
patterns, have a relative difference of only 1.64%. This excellent agreement shows that the approach is not 
sensitive to the orientation of the interface element. 'When the mesh is reasonably fine, the relative difference 
is 2.99% for delamination growth initiation loads and 2.85% for ultimate loads between the two very 
different mesh patterns. This suggests that the front tracing approach is also not sensitive to the mesh pattern 
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once the mesh has a reasonable density. The above observations are also true in the case of a square 
embedded delamination and the approach works well in detecting the delamination front. The proposed 
approach has been shown to be insensitive to the value chosen from the valid range for the interface spring 
stiffness, the orientation of the interface element, and the mesh patterns once a reasonable refinement has 
been achieved. This same case shows that an “enforced front” approach can yield very different results even 
after significant mesh refinement has been achieved. 

Finally, it should be noted that no numerical difficulties have been encountered in the cases 
examined. This indicates that the approach is robust and easy to apply. The proposed front tracing method 
provides the capability of using relatively simple mesh patterns to simulate the moving delamination of an 
arbitrary shape without adapting the mesh as die delamination shape changes. The approach should prove to 
be valuable in future studies of delamination growth. In fact, it will be used in conjunction with the 
development of an approach to modeling through-the-thickness stitching described in the next section. 
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(2) FINITE ELEMENT SIMULATION OF STITCHING AS AN APPROACH TO 
DELAMINATION GROWTH CONTROL 

Introduction 

In laminated composites, the strength is often limited by their low interlaminar fracture toughness, 
which makes them sensitive to delamination. There are many causes of initial delamination and this initial 
delamination can grow into the undamaged area of the composite under certain loading conditions resulting 
in severe strength degradation. In many applications, the delamination alone may be the strength-limiting 
factor, especially if it leads to global instability. Therefore, it is very important to control the delamination 
growth in laminated composites and through-thickness stitching is a very effective way to do this. 

Experimental studies [1-5] show that stitching improves the interlaminar fracture toughness up to 15 
times under certain loading conditions depending on the stitch material and geometry. Due to this reason, 
stitching has already been used by many aircraft and composite structural parts industries. Several stitch 
geometries can be obtained by varying parameters such as stitch materials, fiber/yam diameter, stitch pattern 
and stitch density and the interlaminar fracture toughness depends on the stitch geometry. For example, the 
interlaminar fracture toughness increases in proportion to the stitch density i.e. the number of stitches/unit 
area. For these reasons, it is important to simulate the delamination growth in laminated composites that are 
stitched through their thickness to study the effectiveness of stitches as well as to optimize the stitch 
geometry. 

Lalit K. Jain and Yiu-Wing Mai [6] proposed a micromechanics based model to predict Mode I 
delamination toughness (KIC) for a given stitch parameters using finite difference approach. However, this 
model does not predict the delamination growth length for the applied load or applied opening displacement. 
Instead, the delamination growth value, Aa, has to be assumed and KIC can be found out for that particular 
Aa, for a given stitch parameters such as stitch density, thread material and thread diameter. 

J. Byun et al. [7] tried to simulate Mode I regular and tabbed DCB tests conducted by Guenon [2,3] 
using 2D finite element analysis. Plane stress conditions were assumed. Through-thickness fibers were 
simulated using 3-noded beam elements. Fundamental inputs to the model were the z-axis fiber properties, 
fiber architecture, fiber volume fraction and the Mode I critical strain energy rate (GIC) for the 2D laminate. 
But, this simulation did not model the fiber fracture after elastic stretching. Hence, this analysis cannot be 
considered to be sufficient to model the Mode I delamination growth process. 

Y. Tanzawa et al. [8] modified the model proposed by J. Byun et al. [7] to include the fiber fracture 
and carried out the numerical analysis. As Byun et al. they modeled the DCB specimen with 2D finite 
element mesh. The z-fiber bundle was modeled as a 3-noded rod element with a constant sectional area. This 
model simulates the Mode I delamination growth process reasonably well. However, this model will work 
only for Mode I DCB case with one particular stitch pattern since (a) it does not calculate GII and Gill and 
(b) does not model the specimen width. 

Methodology 

This section of the report proposes an approach to remove the deficiencies presented in the previous 
numerical simulation attempts mentioned above. This proposed approach models the fiber fracture and it is 
able to incorporate any stitch pattern since it models the specimen width using 3D elements. And more 
importantly, the present approach can simulate not only the unidirectional growth of the through-width initial 
delamination in the baseline Mode I case but also the multidirectional growth of an initial delamination of 
arbitrary shape in the mixed-mode case. Moreover, the stitch behavior and the delamination growth process 
have been modeled more realistically by incorporating the unstable crack growth between stitch rows by 
specifying two critical strain energy release rates such as GIC -unstitched and GIC -stitched at corresponding 
positions. Hence, in this present simulation methodology, whenever a stitch row breaks, the crack quickly 
grows and then can be arrested by the next row of stitches as observed in experiments [1-3]. 

The first step in simulating the delamination growth in laminated composites that are stitched 
through their thickness is to study how it is observed in experiments. Previous experimental studies showed 
that the delamination growth process in stitched laminated composites can be either (a) 3-step process as 
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observed by Ogo [1] and Tanzawa et al. [5] or (b) 4-step process as observed by Guenon [2,3] and Dransfield 
et a!. [4], 

Figure 1 shows the 3-step delamination growth process. As shown in (A), a matrix crack passed 
round the thread, and then was arrested at approximately 5 mm past the thread. After loading, the thread was 
stretched and some debonding occurred along the thread as shown in (B). Finally, the thread failed as shown 
in (C) and the stored energy was released to propagate the crack at high speed until the next thread. 



preexiitinj cracks 


debonding 


pull-out miin crack 


Figure 1. 3-steD delamination growth process [11 

Figure 2 presents the 4-step delamination growth process. The first two steps were the same as 
Ogo’s observation. However, in the third step, after debonding, the thread failure occurred inside the matrix 
material and not at the plane of crack. Finally, in the fourth step, the threads completely debonded and pulled 
out from the matrix material. 

The method proposed herein is based on the following assumptions: 

• Fiber fracture is based on Ogo’s [1] 3-step 
delamination growth process as presented in 
Figure i. That is, it has been assumed that trie 
fiber fracture occurs at the plane of delamination 
growth - not inside the matrix material. Hence, 
the friction between the fractured stitch fiber and 
the matrix during pull-out is not modeled. 

• Since the stitch fibers are flexible and they do 
not take up any bending loads, an individual 
stitch has been considered as ID truss element 
and the horizontal connecting fibers as 
shown in Figure 3 are not modeled here. 

• It has been assumed that each individual 
stitch is debonded from the matrix along the 
thread as shown in Figure 1 (B) as well as in 
Figure 2 but is attached to both the top 
surface of the top sublaminate and the 
bottom surface of the bottom sublaminate of 
the composite specimen. This attachment 
simulates the continuity of the horizontal 
connecting fibers that are not explicitly 
modeled here. 




(a) 


(b) 


Figure 3. (a) 3D fiber geometry of orthogonal 
interlocked fabric, (b) Unit cell of the orthogonal 
interlocked fabric composite [2,31 


Interface element 

Xie [9,10] developed an 18-node “interface element” using the UEL option in ABAQUS to calculate 
the strain energy release rate (G) numerically based on 3D-VCCT and to simulate the delamination growth in 
unstitched laminated composites. The advantage of this element is that it can accurately calculate strain 
energy release rate components at an arbitrarily oriented delamination front and the predicted delamination 
front at each increment is not required. In other words, it can predict arbitrary' multi-directional growth with a 
single mesh for the geometry throughout the analysis. The element properties i.e. fundamental inputs for this 
interface element are three spring stiffness values and three critical strain energy release rate (Gi C , Gn C , and 
Gmc) values for unstitched laminated composites (or 2D laminated composites). 

To simulate the delamination growth process in laminated composites that are stitched through their 
thickness, a new “stitch element” is developed using the UEL option in ABAQUS and this new stitch 
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Energy release rate 


element has been combined with the previously developed interface element. Hence, the final finite element 
model for a stitched laminated composite (or 3D laminated composite) specimen, will have two user element 
definitions to simulate the delamination growth along with standard ABAQUS 8-node brick elements that 
model the specimen geometry. 

While defining the element properties of interface elements for this analysis, two different sets of 
critical strain energy release rate (Gi C , Guc, and G ra c) values have to be defined depending on the position of 
those interface elements. For the interface elements that are located in the positions where the specimen has 
z-fiber stitch threads, the strain energy release rate of 3D laminated composites, i.e. Gc-stitched, has to be 
defined. For the other interface elements that are located in unstitched area, the strain energy release rate of 
2D laminated composites, i.e. Gc-unstitched, has to be defined. This will simulate the real life situation, 
which is described as follows. 

During Mode I fracture toughness tests, Ogo observed that the crack growth in unstitched DCB 
specimens was stable. For the stitched specimens, the crack growth was unstable after one row of stitch 
threads was broken, until it got arrested by the next row of stitch threads. That is, the crack growth in the 
matrix region between two rows of stitch threads was unstable. According to Ogo, this is because of the 
different critical energy 
release rate values, Gic, 
between the stitched thread 
area and the unstitched area 
as shown in Figure 4. G IC 
was not constant along the 
beam. Once the threads 
were broken, energy release 
rate G is larger than Gic of 
the unstitched area, i.e. G > 

G IC , which is an unstable 
condition. Then the energy 
was released and the crack 
propagated until the next 
stitching line in an unstable 
manner as shown from 
point A to point B in Figure 
4. If there had not been 
another stitching line, the 
crack would have 
propagated until point C. 

Finally, the loads and boundary conditions are specified to the model and solved using 
ABAQUS/Standard solver. From the output results, the delamination growth area and/or total delamination 
area can be calculated. A general ABAQUS input file to simulate delamination growth processes in . 
laminated composites that are stitched through their thickness will have a format as shown in Figure 5. 



Stitch element 

As observed in experiments, an individual stitch crossing a delamination undergoes elastic stretching 
and eventually fails when its axial stress exceeds the failure strength of the z-fiber yam. A simple “stitch 
element” is developed to represent an individual stitch with this behavior. The formulation of this element is 
based on 2-node truss element. The standard ABAQUS 2-node truss element cannot be used for this purpose 
because it simulates only the elastic stretching and does not simulate fiber fracture. This necessitated the 
formulation and development of the stitch element that models the elastic stretching the same as the standard 
2-node truss element and extends it to have the fiber fracture. 
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Figure S. General ABAQUS input file format to 
simulate delamination growth process in stitched 
laminated composites. 


exceeds the fiber strength. In other words, the element 
stiffness matrix [K] becomes zero when a z > Of where 
CT f is the failure strength of the z-directional fiber 
material. Figure 6 presents the flow chart for the user 
subroutine UEL that simulates the stitch element. The 
subroutine is called by the main input file for every 
load increment. As the applied load increases for each 
increment, the axial stress o z increases. So, whenever 
the axial stress exceeds the failure strength in a 
particular increment, the stiffness matrix becomes zero 
and that element does not bear any load from the next 
increment until the loading process is complete. 


The element is implemented using the 
user element subroutine UEL option in ABAQUS. 
A separate FORTRAN file has been written and is 
interfaced with the main ABAQUS input file. In 
UEL, the element is developed by defining the 
stiffness matrix. The material properties are 
defined as user element properties in main 
ABAQUS input file and these values are passed in 
by ABAQUS whenever it calls the subroutine 
during the solving process. The loads and 
boundary conditions are defined in the main 
ABAQUS input file and this definition accounts 
for the right hand side force vector. The nodal 
displacement vector for this element is found by 
using the ABAQUS/Standard solver. The axial 
strain and axial stress are secondary variables that 
are obtained from the derivatives of primary nodal 
displacement solution in the user subroutine. 

The stitch fiber fracture is simulated by 
setting the element stiffness matrix to zero upon 
failure of the stitch. The maximum stress failure 
criterion has been used for this purpose. The 
fracture occurs when the axial stress of the fiber 



Figure 6. Flow chart for user subroutine UEL that 
simulates stitch element. 
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Mode I delamination growth 

For the baseline case of Mode I delamination growth, Guenon’s [2,3] experimental specimen as shown in 
Figure 7 was modeled in ABAQUS using 8-node incompatible mode (C3D8I) elements. In Figure 7, length 
‘L’ is 152.4 mm, width ‘B’ is 19 mm, thickness ‘2h’ is 4.44 mm and initial delamination ‘ao’ = 25.2 mm. 

Gap elements were defined between the two sublaminates to eliminate the possible layer interpenetration 
after delamination. The stitch pattern as shown in Figure 8 was used for the stitched DCB case. This is the 
stitch pattern that was used by Guenon [3] for 3D laminate experiments. Stitch elements were defined at 
corresponding positions to represent individual stitches. Since Guenon’s stitched specimen is 3D orthogonal 
interlocked woven fabric composite, it has carbon fibers in the z-direction as well. The diameter for the stitch 
element was calculated from the data given in reference [3]. The properties of the specimen and the user 
stitch element are given in Tables 1 and 2. 


BC: u = 0, v = 0, w = 10 mm (unstitched) 



Figure 7. Mode I Doable Cantilever Beam (DCB) specimen. 
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Figure 8. Stitch pattern used in 
stitched DCB model. 


The user interface elements were defined in the initially undelaminated area. The user interface 
element properties that are specified for the unstitched DCB case is presented in Table 3. For the stitched 
DCB case, two G[ C values were specified. G 1C - stitched, which is the fracture toughness value obtained from 
experiments by Guenon [2,3] for a 3D composite laminate was specified for the interface elements that are in 
positions where stitch elements were defined. For the remaining interface elements that exist in the 
unstitched area, Gic - unstitched value was specified. In other words, all the stitch row interface elements 
have Gic - stitched values where as the interface elements between two stitch rows have Gic - unstitched 
values. This was done to simulate the unstable crack growth between stitch rows observed by Ogo [1] as 
discussed before. The Gic - stitched and Gic - unstitched values that were used in the analysis are presented 
in Table 4. The other properties are the same as presented in Table 3 for all the interface elements in the 
stitched DCB model. 

Displacement controlled analysis was performed for both stitched and unstitched cases. A linearly 
increasing vertical displacement ‘w’ was specified for the nodes at the initial delamination end of the 
specimen and the reaction force was requested as an output. The other displacements ‘u’ and ‘v’ in those 
nodes were restrained to simulate the hinge load condition. The vertical displacements ‘w’ for the nodes at 
the undelaminated end surface were restrained. The boundary conditions for both cases are shown in Figure 
7. For both cases, the finite element mesh had 218 elements along x-direction, 8 elements along y- and 1 
each for the sublaminate in the z-direction. Geometric nonlinear analysis was performed with direct control 
time incrementation. 
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Numerical Results 

Figure 9 presents the load- 
displacement curve for the stitched DCB 
model. The delamination growth contour 
plots that show the total delamination 
length ‘a' and the delamination growth 
area ‘AA’ corresponding to the positions 
shown in the load-displacement curve in 
Figure 9 are presented in Figure 10. The 
w'hite area is the initial delamination, the 
cross marks represent broken stitches, and 
the dark gray area represents delamination 
growth. As the applied opening 
displacement increases, the reaction force 
increases steeply up to point 1 since the 
first row of stitches bear the load and 
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Figure 9. Load-displacement curve for the stitched DCB model. 
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Figure 10. Delamination growth contour plots at various load-displacement points - stitched DCB model. 
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elastic stretching of first row of stitches occurs during this period. At point 1, as shown in Figure 10, five 
center stitches in the first row fail and hence there is a sudden drop in reaction force after point 1 . 
Furthermore, at point 2, the remaining two stitches in the first row break and the reaction force drops further. 
Again the reaction force increases up to point 3 and now the second row of stitches carries the load. When 
the applied opening displacement reaches point 3, the five center stitches in the second row fail as in the first 
row. Hence the reaction force drops down to point 4 and at that point, the rest of the stitches in the second 
row fail. This causes another drop in the reaction force after point 4. Similarly, the behavior of the laminate 
and the stitches for other points in the load-displacement curve can be observed from Figure 10. 

The initial delamination in the specimen does not grow until the opening displacement reaches point 
1 1, i.e. 5 = 16.64 mm. From this point onwards, the initial delamination starts growing as shown in Figure 

10. At this instance, it can be observed that the 9 th row of stitches, which are precisely located at the initial 
delamination front, have not failed yet. From Figure 10, it can be seen that these stitches still undergo elastic 
stretching at this opening displacement. Even at position 13, no stitches in that row have failed. Thus the 
method exactly simulates step one of the delamination growth process, i.e. the delamination grows past 
round the stitches and then gets arrested at some distance past those stitches, observed experimentally by 
Ogo [1], Guenon [2,3] and other researchers. Further increases in the opening displacement cause additional 
rows of stitches to fail and the delamination continues to grow in a controlled way. When the opening 
displacement 5 = 40 mm, the total delamination length ‘a’ is 41.3 mm and the delamination growth area ‘AA’ 
is 305.9 mm 2 . 

Figure 1 1 presents the behavior of the edge stitches in each row as the opening displacement 
increases from 0 to 40 m. As the applied opening displacement increases, the axial stress in the first row- 
stitch increases until it breaks. This increase in axial stress is accompanied by the elastic stretching of that 
stitch. When the axial stress exceeds the failure strength ‘of of the stitch material, the stitch fails. In Figure 

1 1, the cross marks again indicate the stitch failure. Figure 12 presents applied opening displacement versus 
total delamination length ‘a’ plotted for the stitched DCB model. It can be observed that the initial 
delamination ‘ao’ does not grow until applied opening displacement reaches the value of 16.64 mm. After 
this point, it quickly grows up to 27.3 mm and is then arrested by the next row of stitches. This shows that 
the simulation behavior conforms to the unstable growth behavior between stitch rows observed by Ogo [1]. 


Opening dspJacwmnt V« Axial tVm on Slitah** - DCB Sttehad modal 



Opening rispiecemert Vs Cr«* length - DCB Stitched model 



Figure 11. Elastic stretching and failure of stitches 
for 6 = 0 to 40 mm. 


Figure 12. Applied opening displacement ‘S’ 
versus total delamination length ‘a’. 


Effects of stitches 

The effects of stitches can be observed by comparing the results of the stitched DCB model with 
those of the unstitched DCB model. The load-displacement curves of both unstitched and stitched DCB 
models for 5 = 0 to 20 mm have been plotted for comparison in Figure 13 where the dashed line represents 
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the unstitched model and the solid line represents the stitched model. For the same applied opening 
displacement, there is a huge rise in the reaction force due to the inclusion of stitches. Hence, it can be 
concluded that the through-thickness stitches greatly improve the z-directional resistance of the specimen. 
Similarly, Figure 14 compares the variation of total delamination value ‘a’ as the opening displacement 
increases from 0 to 20 mm. In the unstitched model, the initial delamination starts growing from the applied 
opening displacement value of 2 mm, whereas in the stitched model, the initial delamination does not grow 
until the opening displacement reaches 16.64 mm. This is due to the improved resistance in the z-direction of 
the specimen provided by the inclusion of the through-thickness stitches. After the initial delamination starts 
growing, the growth is more rapid in the unstitched specimen than in the stitched specimen. When the 
applied opening displacement is at 20 mm, the total delamination in the unstitched specimen is 85.4 mm, 
whereas in the stitched specimen it is 27.3 mm. Hence, it can be concluded that the stitches effectively arrest 
the delamination growth and this simulation methodology works well. 
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Figure 13. Comparison of load-displacement curv es 
between unstitched and stitched DCB models. 
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Figure 14. Comparison of total delamination length 
between unstitched and stitched DCB models. 


Mixed-mode delamination growth 

For the mixed-mode delamination growth. De Moura's [12] carbon/epoxy laminate specimen that 
contains an initial delamination of arbitrary shape created by low-velocity impact was modeled. One quarter 
of the laminate was modeled with C3D8I elements assuming symmetry about both x- and y-axes. The initial 
delamination of arbitrary shape exists at the interface of the top 0° and 90° sublaminates. The schematic 
diagram of the 'A - plate model is presented in Figure 15. The stitch pattern used for the stitched mixed-mode 
case is also shown. Stitch elements were defined at corresponding positions to represent individual stitches. 
The properties of specimen and the user stitch element are given in Tables 5 and 6. 

BC: u = -0.15 mm, w = 0 (unstitched) 



Geometry and applied displacement Stitch Pattern FEM mesh and boundary conditions 


direction shown on V< plate. 

Figure 15. Geometry definition and modeling details for mixed-mode delamination growth problem. 
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Gap elements were defined between the top 0° and 90° sublaminates to eliminate the possibility of 
layer interpenetration after delamination. Also, the user interface elements were defined in the initially 
undelaminated area between these sublaminates. The user interface element properties specified for the 
unstitched case are presented in Table 7. For the stitched DCB case, two different Gic and G n c values were 
specified. Gi C - stitched and Gn C - stitched were specified for the interface elements that are in positions 
where stitch elements were defined. For the remaining interface elements that exist in the unstitched area, G IC 
- unstitched and Gnc - unstitched values were specified. The stitched and unstitched critical G values that 
were used in the analysis are presented in Table 8. All the critical G values were obtained from reference 
[ 12 ]- 


Load-Displacement curve - Stitched 


The finite element mesh and the boundary conditions for both stitched and unstitched cases are 
presented in the third part of Figure 15. The mesh has 20 elements along x-axis, 20 elements along y-axis and 
4 elements through thickness. Displacement controlled analysis was performed, i.e. linearly increasing 
compressive end displacement along x-axis was applied and reaction force was requested as one of the 
outputs. Geometric nonlinear analysis performed with direct control time incrementation. 

Figure 16 presents the load-displacement curve for the stitched mixed-mode model. The 
delamination growth contour plots that show the delamination growth area ‘AA’ corresponding to the 
positions shown in the load-displacement curve in Figure 16 and are presented in Figure 17. As the applied 
compressive end displacement ‘u’ 
increases, the reaction force 
increases almost linearly until 
position 1 in Fi sure 17. At point 1, 
due to the presence of stitches, the 
total laminate starts to buckle 
locally as a plate. Hence, after 
point 1, the slope of the load- 
displacement curve varies. No 
delamination growth occurs until 
point 3, i.e. u = 0.0985 mm. From 
this point onwards, the initial 
delamination starts growing and 
the load-displacement curve 
experiences more nonlinear 
behavior. No stitch has broken 
throughout the analysis. Instead, 
the delamination grew in a mixed- 
mode around the stitches. It has 
been observed that the specimen 
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Figure 16. Load-displacement curve for the stitched mixed-mode model. 
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Figure 17, Delamination growth plots for points 1-10 on the load deflection curve. 
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experiences global wide column buckling from position 9, i.e. when u = 0.203 mm, in the load displacement 
curve. Also, it was inferred that the delamination growth in this stitched mixed-mode model is dominated by 
Mode II, which is in-plane shear. This has been discussed in detail in reference [11], Figure 18 presents the 
applied compressive end displacement ‘u’ versus total (initial plus growth) delamination area ‘A’ plot for the 
mixed-mode stitched model for u = 0 to 0.25 mm. 

The effects of stitches in the mixed-mode 
case can be observed by comparing the results of the 
stitched mixed-mode model with those of the 
unstitched mixed-mode model. The load- 
displacement curves of both unstitched and stitched 
mixed-mode models for u = 0 to 0.15 mm have been 
plotted for comparison in Figure 19. Similarly, 

Figure 20 compares the variation in total 
delamination area ‘A’ as the applied compressive 
end displacement ‘u’ increases from 0 to 0.15 mm 

The plot in Figure 19 shows that the peak 
compressive load increases by a factor of 1 .41 due 
to the inclusion of stitches. Figure 20 shows that in 
the unstitched model, the initial delamination starts 
growing when the applied compressive end 
displacement ‘u’ is at 0.052 mm, whereas in the 
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Figure 18. Applied “u” versus total delamination area. 


stitched model, the initial delamination does not grow until the compressive end displacement value of 
0.0985 mm. After the initial delamination starts growing, the growth is more rapid in the unstitched specimen 
than in the stitched specimen. In the unstitched case, the specimen delaminates completely, which means the 
total delamination area in the !4 - plate is 900 mm2 when the applied end displacement ‘u’ is 0.1045 mm, i.e. 
before ‘u’ reaches its applied maximum value of 0.15 mm. Whereas in the stitched case, the total 
delamination area ‘A’ in this % - plate when ‘u’ is at 0,15 mm is 445.5 mm2. Hence, it can be concluded that 
the stitches effectively arrest the delamination growth and this simulation methodology works efficiently in 
the mixed-mode case as well. 




Figure 19. Stitched and unstitched load-displ curves. Figure 20. Stitched and unstitched delam-displ curves. 


The method can also be used to evaluate the effects of stitching being applied less liberally and 
concentrated in certain geometric patterns. For example, consider the two patterns shown in Figure 21 for 
the quarter plate model. The initial delamination, loading, and boundary conditions are unchanged from the 
previous case. In Figure 22, load deflections plots for these new patterns (#2 and #3) are compared to 
results for the original grid stitching pattern (here identified as #1) and the unstitched plate. It is seen that the 
limited stitching can be used to with good effectiveness to increase the ultimate load. 
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Figure 21. Stitch patterns #2 (left) and #3 (right) applied to the quarter plate modeL 



Figure 22. Load-displacment curves for three stitching patterns compared to no stitching. 


Conclusions 

The purpose of this research was to simulate the delamination growth process in laminated composites that 
are stitched through their thickness using finite element analysis so that this simulation can complement 
experimental techniques used to study the effectiveness of stitches as well as to optimize the stitch geometry 
by varying the parameters such as stitch materials, fiber/yam diameter, stitch pattern and stitch density for 
any given loading case. ABAQU S/Standard was used for this purpose and a simple “stitch element” was 
developed using the UEL user element subroutine option. This “stitch element” was combined with a 
previously developed 3D-VCCT based algorithm that can trace the moving front of a delamination of 
arbitrary shape. The proposed methodology was successfully applied for both the baseline case of Mode I 
delamination growth and the mixed-mode delamination growth. The results showed that this methodology 
works effectively and simulates the behaviors observed in experiments. Due to the simplicity of this 
methodology, i.e. modeling an individual stitch and placing it in the laminate at corresponding positions 
depending on the stitch pattern, various stitch materials, patterns, and densities can also be easily evaluated to 
optimize the stitch geometry depending on the applications and/or boundary conditions. 
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Table 1. Carbon/Epoxy laminate properties. 


Table 5. Carbon/Epoxy laminate properties. 


Properties 

Values 

E„ 

21.07 GPa 

E 22 

20.76 GPa 

E 33 

8.27 GPa 

V12 

0.0357 

V13 

0.0426 

V23 

0.0426 

G« 

2.5 GPa 

G,3 

2.0 GPa 

G 23 

2.0 GPa 


Table 2. Stitch element properties. 


Properties 

Values 

Cross section area, A 

0.151 mm 2 

Young’s modulus, E 

228.0 GPa 

Poisson’s ratio, v 

0.21 

Failure strength, <jf 

3.45 GPa 


Table 3. Interface element properties. 


Properties 

Values 

K, 

1 .0 E+07 N/mm 

k 2 

1 .0 E+07 N/mm 

K, 

1.0 E+07 N/mm 

Gic 

0.307 N-mm/mm 2 

Ghc 

0.632 N-mm/mm 2 

Ginc 

0.817 N-mm/mm 2 


Table 4. Interface element properties. 


Properties 

Values 

GIC - stitched 

3.850 N-mm/mm2 

GIC - unstitched 

0.307 N-mm/mm2 


Properties 

Values 

Ell 

109.34 GPa 

E22 

8.82 GPa 

E33 

8.82 GPa 

vl2 

0.342 

vl3 

0.342 

v23 

0.52 

G12 

4.32 GPa 

G13 

4.32 GPa 

G23 

3.20 GPa 


Table 6. Stitch element properties. 


Properties 

Values 

Cross section area, A 

0.151 mm 2 

Young’s modulus, E 

235.0 GPa 

Poisson’s ratio, v 

0.21 

Failure strength, a f 

3.73 GPa 


Table 7. Interface element properties. 


Properties 

Values 

K, 

1.0 E+07 N/mm 

k 2 

1.0 E+07 N/mm 

K 3 

1.0 E+07 N/mm 

Gic 

0.306 N-mm/mm 2 

Gnc 

0.632 N-mm/mm 2 

Gmc 

0.817 N-mm/mm 2 


Table 8. Interface element properties. 


Properties 

Values 

Gic - stitched 

3.850 N-mm/mm 2 

Gnc - stitched 

6.320 N-mm/mm 2 

Gic - unstitched 

0.306 N-mm/mm 2 

Gnc - unstitched 

0.632 N-mm/mm 2 
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(2) Damage Tolerance of Mixed-Tow Hybrid, Composite Laminates 

John M. Kennedy, Professor of Mechanical Engineering, Clemson Univ. 

Josh Stewart, Graduate, Mechanical Engineering, Clemson Univ. 


Summary 

The objective of this research is to investigate damage tolerance and durability 
concepts for mixed-tow hybrid (MTH) composites as well as buffer strip panels. An 
experimental study on MTH composites was completed in the first two years of this 
project. The data from the experimental program provide the basis for developing a 
model that can predict the failure behavior of MTH composites. Presently, a two- 
dimensional finite element code is being developed to predict damage accumulation and 
progression in a notched buffer strip panel subject to tensile loading. The model must be 
capable of identifying stress levels in each ply of the laminate as well as predict ply level 
failures near stress concentrations. The models, once completed, will account for 
material inhomogeneities, anisotropy, and to some extent characteristic damage. 

Another objective of this research is to experimentally examine the combined 
effect of buffer strips or MTH configurations with through-the-thickness (TTT) 
reinforcement. Studies in which both techniques have been used together have been 
limited. This project involves carbon fiber material systems with both techniques 
employed. 

Section I: Narrative 

The use of mixed-tow or intraply hybrids as a damage tolerance concept has been 
investigated at Clemson University. Research is being conducted under two subsets as 
described below: 

(1) Finite Element Analysis of Buffer Strip Panels 

An experimental investigation of damage mechanisms in mixed-tow hybrid, 
center-cracked tension panels was conducted. Five mixed-tow combinations were 
studied. All mixed-tow plies in the panels were made with uniweave cloth containing 
periodic placement of glass and carbon tows. A quasi-isotropic and 0°-dominated 
laminate were evaluated. Results of radiographs and COD measurements indicated that 
the introduction of glass fibers into carbon fiber panels improved the fracture resistance 
of those panels. A decreasing carbon tow/glass tow ratio cause the predicted critical- 
damage zone size to become larger; therefore, the critical strain intensity factor increased, 
and the panels were less notch sensitive. 

Damage tolerance and durability criteria for metals allow damaged structures to 
be in service until the damage reaches a specified size determined by the critical damage 
size. Only when a similar durability criterion for composites is established will it be 
feasible for composites to be the primary structural material in aircraft. The bridge to 
achieving that standard lies in the incorporation of novel material design concepts other 
than using excessive material, which increases cost. 

Mixed-tow hybridization (often referred to as intraply hybridization) has 
been shown to improve damage tolerance characteristics over all-carbon systems in 
tension-loaded composite laminates [1-5]. In [1], it was confirmed that notch 
sensitivity of a material might be reduced by mixed-tow hybridization with higher 
failing strain fibers (i.e., S-Glass), as had been previously observed in buffer strip 
panels [6,7], In addition to the fiber with higher failing strain, this reduction in 
notch sensitivity was attributed to the extensive damage that occurred in the notch- 
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Figure 1. Schematic of mixed-tow hybrid laminate. 


tip region prior to final fracture. The introduction of the higher strain and lower 
modulus fibers into the carbon material system induced more pre-catastrophic 
damage. A mixed-tow hybrid involves a repeating pattern of a hybridizing tow or 
tows in some or all plies of a laminate as shown in Figure 1 . The hybridizing tows 
can be readily introduced into the composite fabric in a weaving process. 

In [2], S2-glass was used as the hybridizing material since the buffer strip 
studies found S2-glass buffer-strip material gave the most damage tolerance 
improvement over all-carbon panels [6]. Further, including hybridizing tows in all 
plies resulted in a higher strength and fracture toughness than including tows in 
only the primary load-bearing plies for panels with central slits. 

The objective of this paper is to examine the effect introducing hybridizing high 
failing strain glass tows has on damage mechanisms and consequently, fracture 
properties of center-cracked tension panels. Two laminate orientations and several 
mixed-tow combinations were investigated which were not investigated in [1-5]. 
Mechanisms that result in improved fracture toughness of mixed-tow hybrid 
composites were identified. The development of subcritical damage at the crack tip 
was studied using x-radiography. Panel strengths are compared to predictions from 
a linear elastic fracture mechanics (LEFM) model. 


EXPERIMENTAL PROCEDURES 
Materials and Specimens 

The specimens were made with IM7-6K carbon manufactured by Hexcel and 
449 S-2 Glass fiber roving manufactured by Owens Coming. The fibers were 
woven into a uniweave composite fabric using a 220-denier polyester fill yam 
manufactured by Trevira of Hoechst Celanese. Cloth was woven with five mixed- 
tow combinations: all-carbon, all-glass, and three hybrid system ratios with carbon 
tow to glass tow ratios of 10:2, 6:2, and 2:2. Figure 2 shows pictures of the 
uniweave cloth with the five mixed-tow combinations. The tow size was such that 
there were about seven tows per centimeter. The resin system used in the material 
was 3501-6 manufactured by Hexcel. Ten sheets were manufactured at NASA 
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(a) All-carbon 


(e) All-glass 


(b) 10:2 (c) 6:2 (d) 2:2 

Figure 2. Composite fiber mixed-tow combinations. 


Langley Research Center with a constant mass fraction of 32%. Thus, the fiber 
volume fraction, Vf, varied due to the different densities of the carbon and glass. 
All results have been normalized to a common fiber volume fraction. 

Two laminate orientations were manufactured using uniweave cloth for every 
ply. One laminate was a quasi-isotropic layup, [+45/0/-45/90]s, representative of 
that used in fuselage structure, while the other was a 0°-dominated layup, [+45/- 
45 / 02 / 90 / 02 /- 45 /+ 45 ]t, representative of that used in wingskin structure. Test panels 
were cut from the sheets to the sizes shown in the test matrix in Table 1. Two sizes 
of fracture panels were tested with the notch lengths indicated in the table. Two 
replicate tests were conducted for each configuration. 

Test Procedures and Equipment 

All tests were conducted at room temperature in a 250-kN servo-hydraulic 
testing machine. Panels were loaded to failure at a loading rate of 8.9 kN/min. 
Load, crack -opening displacement (COD), and strain (selected samples) were 
recorded using a computer-based data acquisition system. COD data was measured 
with a clip-gage. A steel fixture was clamped to the small panels to prevent out-of- 
plane displacement at the center of the notch that minimized Mode III loading at the 
notch tips. For the large panels, an aluminum fixture that spanned the width of the 
specimen was used to prevent out-of-plane displacement. 


Representative results are presented since the large number of tests conducted 
prevents the presentation of all test results in this paper; however, complete results 
may be found in [9] including those of unnotched coupons. 

Fracture Tests Evaluation 

Figure 3 shows the basic COD behavior of the mixed-tow hybrid panels. The 


TABLE I. SPECIMEN GEOMETRIES 



2a/W 


0.063 


0.625 

0.313 

0.625 


304.8 


304.8 

609.6 

609.6 

Specimen width (mm): ! 

101.6 

101.6 

101.6 

203.2 

203.2 

Notch length, 2a (mm): 

6.4 

25.4 

63.5 

63.5 

127 




















glass in the panels caused arrest-crack behavior represented by the jumps in the 
curves, much like buffer strip panels. The panels with longer notches generally had 
more COD and consequently, more precatastrophic damage in the notch-tip region. 
Also, the 203.2-mm wide panels with a 63.5-mm notch underwent more COD than 
the 101.6-mm wide panels with the same notch length. There was also more pre- 
catastrophic damage in these panels. 

Typical radiographic and COD results are shown in Figure 4 for a mixed-tow 
hybrid panel. The black bar in the middle of the radiographs is the out-of-plane 
displacement constraint mentioned previously. The damage growth indicated in the 
radiographs was quantified by sudden increases in COD. 

Figure 4 shows the significance of the initial notch-tip location relative to the 
glass tows and the impact it has on damage modes and the damage-zone size. The 
first radiograph, taken at 58% of failing stress, shows that little damage developed 
at the notch tips prior to a significant jump in COD. The left notch tip was adjacent 
to glass tows in one of the primary load-bearing (0°) plies. At higher stress, some 
significant matrix cracking in the 90° plies is evident with some delamination close 
to the notch tip. There is also minimal cracking in the 45° plies, while there is a 
large axial split parallel to the 0° plies. Conversely, for the right notch tip that was 
located some distance away from the 0° glass tows, there were no axial splits. The 
damage developed as matrix cracks in the 45° plies and delamination. The damage 
did not extend outboard of the 0° glass tows. 
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Figure 3. Stress vs. COD for 6:2 ratio, [+45/-45/0 ; /90/0 ; /-45/+45] T panels. 
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Figure 4. Stress vs. COD and radiographs (6:2 ratio, quasi-isotropic, 4-inch wide, 1-inch notch). 

Figure 5 shows that the damage-modes and size at both notch tips in the 2:2 
ratio panel were very similar. For these panels the initial location of the notch tip 
was less significant and the damage in the panel was symmetric, similar to the 
behavior observed in the all-carbon and all-glass systems. The close spacing and 
the fact that the glass tow architecture was not aligned in plies of the same angle 
resulted in approximately the same amount of glass in both notch-tip regions. 

The presence of glass in the notch-tip region generally increases the amount of 
pre-catastrophic damage. Figure 6 shows evidence of this trend with mixed-tow 
combinations of 10:2 and 6:2 at 92% of failing stress. The damage-zone for the 6:2 
ratio is larger at both notch tips. There are large axial splits along the 0° and 45° 
plies accompanied by extensive delamination. Pre-catastrophic damage in the 
notch-tip region redistributes the stress away from the notch tip; thus the panel 
becomes less notch sensitive. The radiographs of the 6:2 ratio panel suggests that 
damage development was affected by glass tows at the notch tip and the next set of 
tows outboard from the notch tips. This is not evident in the 10:2 ratio panel. The 
reason for the large amount of damage in the 6:2 panel is that the failing strain of - 
glass fibers is two to three times higher than that of the carbon fibers and the matrix 




wide, with 25.4-mm notch. o.=444.7 MPa 


I 
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material, which have about the same failing strain. Therefore, the matrix fails while 
the glass fibers stay intact, whereas in the all-carbon panel, the matrix and the fibers 
fail at about the same strain levels. This explains the general trend in the panels for 
which pre-catastrophic damage increased as the carbon tow/glass tow ratio 
decreased. 

The increased strain that occurs in the glass tows during a significant damage 
event was evident from the failing strains of the mixed-tow hybrid panels. Figure 7 
shows that the failing strains of the panels increased with a smaller spacing of the 
glass tows, which also meant there was more glass in the panels. For each notch 
length, the failing strain was generally highest for 2:2 ratio panels and lowest for 
10:2 ratio panels. 


Fracture toughness was determined for each panel using a LEFM model 
developed for composites [8]. The stress intensity factor equation for a centrally 
cracked sheet of infinite width, K a =o^Ka , was the basis for these calculations. 
A critical damage-zone size p c was incorporated into this equation such that the 



(a) 10:2 (b) 6:2 

Figure 6. Quasi-isotropic 203.2-mm wide panels with 63.5-mm notches. 
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Figure 7. Failing strain results for [+45/0/-45/90] s 
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critical stress intensity factor, Koq became 


K oQ =o c Jn(a + p c ) (1) 

where o c is the far-field stress and a is the initial slit half-length. Then p c , assuming 
Pc> 0, reduced the predicted failing stress by assuming a longer crack than the initial 
machined slit. However, a complex damage-zone did exist in the notch-tip region, 
and using p c to account for the effect of the damage-zone on failing stress was 
adequately effective. 

The critical damage-zone p c can be determined for a notch length of zero where 
a is equal to the material ultimate tensile strength, F m . Thus, Koq was used to 
obtain p c as follows. 


Pc = 


<%> 


n 


( 2 ) 


To determine the critical stress intensity factor for the panels tested, (2) was 
substituted into (1) to yield 


*«=■ 
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7Easec 
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sec 


Jta 

W 


( 3 ) 


Included in this equation is a widely used finite-width correction factor 
yjsec(7ta/W) where ITis panel width. 

The critical stress intensity factor was averaged for all the panels of a mixed-- 
tow combination and laminate orientation. This provided a universal stress 
intensity factor for each mixed-tow combination and laminate orientation. Fracture 
is also characterized in terms of a strain intensity factor. The critical strain intensity 
factor, Keq, was determined by dividing Kcq by the material’s elastic modulus. The 
fracture properties of each mixed-tow combination and laminate orientation are 
presented in Table 2. 

The results in Table 2 show that there was no identifiable trend relating the 
carbon/glass tow ratio to the critical stress intensity factor. In terms of the stress 
intensity factor, the all-glass laminates exhibited the best fracture toughness 
followed by the 2:2 ratio. For the 0°-dominated lay-up, the all-carbon laminates 
had the third best K^q followed by the 10:2 ratio and the 6:2 ratio. For the quasi- 
isotropic lay-up, the 10:2 ratio had the third best Koq followed by the all-carbon and 
the 6:2 ratio. However, the strain intensity factor of the all-carbon material is the 
lowest and that of the all-glass is the highest. The decreasing elastic modulus with 
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TABLE 2. AVERAGE FRACTURE PROPERTIES. 


Laminate orientation 

Mixed-Tow 

Combination 

Average Koq 
(MP aVmm) 

Average Kjq 

(Vmm) 

[+45/-45/0 2 /90/0 2 /-45/+45] t 

All-carbon 

4451 

0.05519 

10:2 

4098 

0.05715 

6:2 

3515 

0.05635 

2:2 

4466 

0.07930 

All-glass 

5810 

0.17585 

[+45/0/-45/90]s 

All-carbon 

3160 

0.06209 

10:2 

3224 

0.06789 

6:2 

2879 

0.07081 

2:2 

3625 

0.09762 

All-glass 

3680 

0.17529 


increasing percent glass in the laminates resulted in higher strain intensity factors 
for the mixed-tow combinations. 

Although the all-glass laminates had the best fracture toughness, one of the 
primary benefits of carbon fiber composites must be kept in mind, that being the 
lightweight characteristic. This is an important consideration in the design of 
structural components for aircraft. Therefore, the density of each mixed-tow 
combination and laminate orientation was experimentally determined [9]. Realizing 
that the S-glass has a higher density than the carbon suggests that the best fracture 
toughness in terms of both stress and strain intensity factors may not make it the 
best material choice. Figure 8 displays bar graphs of the fracture toughness of each 
mixed-tow combination and laminate orientation. In this figure, the magnitude of 
each stress or strain intensity factor has been normalized by the density of its 
corresponding material’s density. In addition, the data in each plot was normalized 
to the largest magnitude (the all-glass, 0°-dominated results) so that the fracture 
toughness parameters of the mixed-tow combination and laminate orientation that 
exhibited the highest magnitude were equal to unity. 

Normalizing the fracture toughness parameters by the densities decreased the 
difference in magnitude among the laminates. The stress intensity factor in the 
quasi-isotropic laminates is now highest in the 2:2 ratio mixed-tow combination, 
with that of the all-glass system being next to last instead of the highest as it was 
before density-normalizing. The 0°-dominated laminate orientation still had the 
highest stress intensity factor in the all-glass system followed by the all-carbon 
system instead of the 2:2 ratio as before. 

A decreased difference in critical stress intensity factors was noted in [8] in the 
study of carbon laminates. It was found that K^q was the same for all lay-ups with 
E x = E y . In the current study, Kqq for the quasi-isotropic laminate (Table 2) varied 
by less than 10% for all mixed-tow combinations except the 6:2 ratio, which was 
approximately 13% below the mean. The variation among the mixed-tow 
combinations was even less when Kqq was normalized by the density. There was 
significantly more variation in the 0°-dominated laminates; difference between the 
highest stress intensity factor (all-glass) and the lowest (6:2 ratio) was 40%. 
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Figure 8. Fracture toughness parameters normalized by the density of the corresponding mixed- 
tow combination and laminate orientation with the highest magnitudes being unity. 


The 0°-dominated laminate orientation generally had a lower strain intensity 
factor than the quasi-isotropic laminate orientation, indicating the 0°-dominated 
laminates are more notch sensitive than the quasi-isotropic laminates. This 
correlates \v ith conclusions in [8] that an increased proportion of 0° plies results in 
an increased notch sensitivity of the laminate. 

Using the average values from Table 2, the predictions from the model were 
compared to the experimental results. Solving (3) for cr c gives 
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for finite-width panels. Plots of far-field failing stress versus notch length (2a) 
were generated from this equation for each mixed-tow combination and laminate 
orientation. Plots were also generated from 



for infinite plates. Figure 9 shows stress versus crack length predictions for the 0°- 
dominated, 6:2 ratio laminate. These results are typical of those for all the mixed- 
tow combinations. Predictions from (3), for two finite width plates, a 101.6-mm 
wide plate and a 203.2-mm wide plate are shown. Also included on the plot are 
data actual failing stresses of the panels. Figure 9 also contains two additional 
curves. These curves (dashed lines) represent the net section failing stress for the 
101.6-mm and 203.2-mm wide panels. In all cases, the net section failing stress 
was higher than the failing stress of the panels, thus confirming the stress reduction 
caused by the stress concentration at the notch tips. The LEFM model predicted a 
higher failing stress than the net section model for large 2a/W. This phenomenon is 
a result of the finite-width correction factor in the LEFM model. 

The model predicts the failing stresses of the 203.2-mm wide panels well, but is 
not as accurate for the 101.6-mm wide panels. This inaccuracy is due to the stress 
intensity factors being generally much larger for the 203.2-mm wide samples. 
Since the stress intensity factor for the overall mixed-tow combination and laminate 
orientation is an average value, it was skewed upward a little by the values of the 
203.2-mm wide panels. The cause of the stress intensity factor elevation for the 
larger panels stems from the large amount of crack-tip damage, discussed earlier. 

This behavior is shown graphically in Figure 10, which displays the critical 
stress intensity factor for each of the 10:2-ratio mixed-tow panels. The figure 
shows the trend of panels with longer notches having higher stress intensity factors, 
while larger panels with longer notches had even higher stress intensity factors. A 
variation from the latter trend is evident for the 0°-dominated panels with a 63.5- 
mm notch where wider 203.2-mm wide panels had lower stress intensity factors 
than the 101.6-mm wide panels. This was observed in most of the mixed-tow 
combinations at the 63.5-mm notch length for the 0°-dominated laminate 
orientations. Also, there is a significant amount of variation among the tests, even 
among those of the same laminate orientation and notch length. This variation may 
not be surprising as the location of the notch tip relative to the 0° glass tows 
strongly affected pre-catastrophic damage. In light of these results, Figure 9 
indicates the fracture strengths of the 203.2-mm wide laminates are more accurately 
predicted than those of the 101.6-mm wide laminates. Walker, et al discussed in [2] 
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Figure 9. Failing stress versus crack length for 6:2 ratio, 0°-dominated panels. 
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Figure 10. Stress-intensity factor versus 2a for 10:2 ratio laminates. 


the likelihood that “methods based on the classical square-root stress-field 
singularity are not capable of predicting strengths for a large range of notch sizes, 
and that reduced-singularity methods may be more accurate.” However, the LEFM- 
based predictions in the current study are within reasonable accuracy 


The effect on fracture toughness of weaving various ratios of hybridizing S- 
glass tows into a carbon fiber material system was studied. Damage mechanisms 
were examined through study of crack opening displacement and x-radiography, 
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while fracture-toughness parameters were determined using LEFM. The following 
general conclusions were reached from the results of this work. 

• The critical strain intensity factor increased with the glass/carbon ratio in 
laminates; thus causing the laminates to become less notch sensitive. 

• No trend was identifiable that related the carbon/glass ratio to the critical stress 
intensity factor. 

• Panels with larger notches experience more pre-catastrophic damage than panels 
with shorter notches. 

• Pre-catastrophic damage in the notch-tip region reduced notch sensitivity and 
increased the strain intensity factor. 

• Alignment of the glass tow architecture among 0° plies along with initial notch- 
tip location relative to the nearest glass tows influenced non-critical damage. 

• Mixed-tow panel crack-arrest behavior is like that of buffer strip panels. 

• Using the critical stress intensity factor enabled the prediction of the failing 
stress in the laminates within reasonable accuracy. Generally, strength 
predictions for panels with short notches were unconservative. 
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(2) Combined Effect of Buffer Strips or MTH Configurations with TTT 
Reinforcement 

The objective of this research was to examine the effect of R-ratio on failure 
modes for through-the-thickness (TTT)-reinforced composite materials. It was 
anticipated that as the R-ration moved from one that was compression-dominated to one 
that was tension-dominated, the failure mode would at some point transition respectively. 
Tests were run in the range from predominaently tension to predominantly compression, 
and all samples failed with some form of longitudinal splitting around the open-hole, 
transverse cracks growing from the open-hole region to the outer boundaries, and 
delamination in the region local to the hole. Photomicrographs from a compression- 
dominated test showed local compression failures in plies and delamination. 

As expected, the IT "1 -reinforcement decreased the intensity of the delaminations 
that occurred. Damage growth was more dependent on stress amplitude than on the 
maximum tensile or the maximum compressive stress. For instance, the fully-reversed 
tests, with a larger stress amplitude, always failed much earlier than the predominantly 
tension-dominated test, with a smaller stress amplitude, even though the maximum 
tensile stress was the same. Similar behvaior was observed with the compression- 
dominated tests; the fully-reversed tests failed much earlier than the predominantly 
compression-dominated test, even though the maximum compressive stress was the same. 
A larger stress amplitude caused localization of damage mech anis ms around the open- 
hole area, whereas smaller stress amplitudes caused a more dispersed region of 
longitudinal and transverse splitting, sometimes as far as two hole diameters away. 

Lastly, secant modulus degradation proved to be a good indicator for estimating 
fatigue life. Although all samples exhibited a gradual decrease in secant modulus, failure 
was imminent when a large and sudden decrease in secant modulus was observed. 

Section II: Measurables 

Publications: 

Lagrange BA and Kennedy JM. "Damage Mechanisms in Mixed-Tow Hybrid, Center- 
Notched Tension Panels," No. 91 in Proceedings of the American Society of Composites, 
16th Conference, Blacksburg, VA (2001). 

Osborn, C.J. and Kennedy, J.M. "Progressive Damage in Mixed-Tow Hybrid Center- 
Notched Composite Tension Panels," No. Proceedings of the American Society for 
Composites 17th Technical Conference, West Lafayette, IN (2002). 


Section IV: Participants: 

Graduate Student: 

Josh Stewart (W, M) 
Faculty: 

John M. Kennedy 


Section V: Collaborations: 

NASA Centers: 

NASA Langley Research Center, Hampton, VA 
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(3) A Three Dimensional Fatigue Crack Closure Model of Surface Cracks 

Paul Joseph, Associate Professor of Mechanical Engineering, Clemson Univ. 

Magda Capitaneanu, M.S., Graduate, Engineering Mechanics 
Adrian Loghin, Ph.D. Graduate, Engineering Mechanics 

Summary 

The original objective of this research was to develop a model to predict fatigue crack growth in the 
fuselage of metallic aircraft structures. A realistic model must include crack closure effects induced by 
plasticity and possibly other nonlinear issues that occur with fracture in metals. Due to the nonlinear, three- 
dimensional nature of the problem, it is highly desirable to use a simplified model with the ability to account 
for these nonlinear effects at every time step. The line spring model was chosen as the basis for this model. 

In building this model, several aspects of linear and nonlinear fracture have been investigated. 

In order to account for a loss of K-dominance, a method was developed to accurately determine the 
coefficients of the higher order terms in the expansion of the stress field near the tip of a crack. For the case 
of an edge crack in a strip, which is the basic geometry required to implement the line spring model, the first 
four coefficients have been determined, where the first is the stress intensity factor, the second is the so- 
called T-stress, etc. A journal article describing the method and results is currently in preparation. The 
symmetric case has been completed, but due to our recent interest in mixed mode fracture, which is discussed 
below, the antisymmetric case will also be included in this article. A female M.S. student, Magda 
Capitaneanu, completed her thesis on this subject and graduated in May of 2000. The symmetric results will 
be implemented into the model to help approximate constraint variation along the crack front. In order to 
make use of the model to predict surface crack growth; it is necessary to extend the capability of the line 
spring model to characterize crack tip conditions near the endpoints of the crack, as well as in the center of 
the crack. Another accomplishment of this research has been the determination of the precise nature of the 
stresses and crack opening displacement predicted by the line spring model at the endpoints. It has been 
determined that it is necessary to include a second singular term in the numerical solution procedure. The 
current method that has been used for nearly thirty years is to use only one singular term. Given this result, it 
will now be possible to devise a method to approximate crack tip conditions all along the curved crack front. 

In the last three years we have extended our work with higher order terms from linear material 
behavior to nonlinear material behavior. By far the most interesting work has been with higher order/mixed 
mode analysis of nonlinear material behavior. Adrian Loghin did his doctoral work on this problem, and 
graduated in August of 2001. We have discovered that multiple asymptotic solutions exist. The important 
point is that the stress-based, J-integral approach fails near mode I, and a CTOD approach might be better. 

Section I: Narrative 

Description of Research and Accomplishments - this section is divided into two parts: 
nonlinear mixed mode fracture and higher order linear elastic terms. 

1. Nonlinear mixed mode fracture 

The primary accomplishment of this research has been on understanding the effects of arbitrary 
loading on the asymptotic stress and strain fields in nonlinear fracture for a stationary crack tip. In this work 
the material behavior is modeled as pure power law hardening within the theory of J 2 -deformation plasticity. 
Of particular interest is the role of two-parameter higher order solutions in “mixed mode” fracture, where 
each of the first two terms in the expansion of the displacement fields is variable-separable with real, but 
unequal eigenvalues. Such an asymptotic form is given below. 

= K, (i7r,) X| u|(0) + K 2 (r/r,) x = u?(0) + ... , \ < X 2 


= a A° (r / r, )“' +1 u, 1 (0) + a A°~' A 2 (r / 


r ) n s l+ l + ( s 2 -» 1 )-2( 0 ) + , 


(3-D 
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e s (r,0) 


= K, (r / r t ) - ^ 1-Xi ^ Ey (0) + K 2 (r/r,)- (l ' Xj) e?(0) + . . . 


(3.2) 


= aA“ (r / r, )“' e y (0) + aA” _1 A 2 (r / r, )“' + ^ 2 ~ s ‘ (0) + . . . , 

aii (r,0) -I I -l(i-x.) , -I I-i -To-*., , 

— - = a "K, n (r/r,) " a*(0) + a a K? K 2 (r/r,) L ^ " ;j oJ(0) + ... 


= A,(r/ii) s ‘ a y (0) + A 2 (r/r,) S2 a y (0) + . . . , s, < s 2 , i, j = r, 0. 


(3.3) 


where Ki and K 2 (or Ai and A 2 ) are arbitrary constants, Xi and X 2 (or Si and s 2 ) are the eigenvalues with 0 < Xi < X. 2 < 1 
(-1 < Si < s 2 < 0) that govern die singular behavior of stresses and strains at the tip of the crack (or wedge), and the 
barred functions are the eigenfunctions which define the angular dependence. 

A displacement-based finite element method coupled with singular value decomposition is used to 
solve the nonlinear higher order eigenvalue problem for cracks and wedges. The versatility of this method to 
account for arbitrary loading and geometry is used to investigate homogeneous crack and interface crack 
cases by studying the limit as a wedge subjected to arbitrary loading becomes a crack. Two such limits are 
considered. The first when the total wedge angle is smaller than the crack case, and the second when it is 
larger. By using this approach, it has been determined that in homogeneous plane strain and plane stress 
fracture, and in plane strain fracture of an interface crack between a hardening material and a rigid substrate, 
there are two distinct asymptotic solutions. These solutions are referred to as mode I dominant and mode II 
dominant. In the linear elastic limit the solutions are identical, however, as the hardening exponent deviates 
from the linear elastic case, the asymptotic structures are different. The forms of the asymptotic solutions are 
presented in the table below. 

Table 3.1. Asymptotic behavior of mixed mode loading cases in nonlinear fracture 
based on Equations (3.1 - 3.3), as determined by the current research 


(Log 

lin and Joseph, 2001). 




Plane strain 
homogeneous 
fracture 

Plane strain 
interface fracture 

Plane stress 
homogeneous 
fracture 

Mode I 
dominant 
loading 

Higher order, s, <s 2 

Higher order, s, < s 2 
Shaima and Aravas (1993) 

Higher order, s 2 <s, 

Mode II 
dominant 
loading 

Mixed mode, s, = s 2 
Shih (1974) 

Oscillatory 

Higher order, s t <s 2 


As seen from the table, each of the three fracture cases has one higher order solution where Si < s 2 . In both 
plane strain cases this higher order solution is mode I dominant, while in plane stress the solution is mode II 
dominant. The mode I dominant interface crack higher order solution was first presented by Sharma and 
Aravas (1993). The “mode II dominant” solution for homogeneous plane strain fracture, first presented in a 
classic study by Shih (1974), is a mixed mode solution governed by a single real eigenvalue. Contrary to this 
mixed mode solution, there is strong evidence that the mode II dominant plane strain interface crack solution 
is oscillatory, i.e., non-separable with a form similar to that of the compressible linear elastic case that results 
from a complex eigenvalue. The mode I dominant plane stress is solution is the most unusual, as it is in a 
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higher order form, but the antisymmetric portion is more singular that the symmetric portion. This unusual 
result is discussed in more detail later in the report. The plane strain case has been published (Loghin and 
Joseph (2003)) and the plane stress case will be submitted for publication soon. The key asymptotic results, 
which reveal these multiple asymptotic solutions, are given in the next three figures. The important point is 
that the J-integral term does not account for the mix of the loading in some cases. These cases include: Mode 
I dominant for plane strain (Figure 3.1), all loading in plane stress (Figure 3.2) and mode I dominant loading 
for the interface crack (Figure 3.3). 



330 340 350 360 370 380 390 

a (degrees) 

Figure 3.1. The first two eigenvalues for a linear elastic (n = 1, v = 0.5) 
and a nonlinear elastic (n=3) free-free wedge. The nonlinear case shows 
multiple asymptotic solutions (see Eqn. 3.3) for the eigenvlaues for the case 
of a crack. 
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Figure 3.2. Same as Figure 3. 1 for plane stress. 
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Figure 7a. 

n=l, linear clastic. 

elampcd-frcc wedge 


More symmetric branch — - 

Mode 1 dominant 


Mere antisymmetric branch 


Mode II dominant 


90 200 


More symmetric branch - 



Mode I dominant 


I Figure 7b. 

' n=1.8, nonlinear elastic, 

plane strain 

Eigenfunctions 
I are identical — . 


Mode n dominant Sl 


a (degrees) 

Figure 3.3. The first two eigenvalues for a linear elastic (n = 1, v = 0.5) 
and a nonlinear elastic (n=3) clamped-ffee wedge for plane strain loading 
conditions. A variable-separable solution cannot be obtained in the lower 
figure between a = 180° and the angle at which the double root occurs. 
Again the figure shows multiple solutions for nonlinear material behavior 
for an interface crack (a = 180° ). 










The results in Figures 3.1 - 3.3 are asymptotic, and should be verified by full-field results. For the cases of 
plane strain and plane stress, we have confirmed the asymptotic results using full field finite element 
analysis. The asymptotic results for all cases predict that a switch in asymptotic structure takes place 
somewhere between pure mode I and pure mode n. For the case of plane strain this is demonstrated in 
Figures 3.4 and 3.5. In Figure 3.4 there is proof that near mode I, somewhere between a loading 
corresponding to 0.8 < M* < 0.98, where NT is the loading parameter defined by Shih (1974), 


M e = — tan" 

7Z 


K, 


K, 


(3.4) 


the asymptotic solution switches from mixed mode where symmetric and antisymmetric parts of the stress 
field have the HRR eigenvalue, to a higher order form where the antisymmetric part of the stress field has a 
weaker eigenvalue. The complicated way in which this occurs is demonstrated in Figure 3.5. For more 
details, see the paper by Loghin and Joseph (2003). 
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Figure 3.4. Comparison of the asymptotic solutions with the symmetric 
and antisymmetric parts of the full-field finite element results for n = 3 at 
0 = 90° for different loading conditions. Figure 3.4a: M e = 0.98, Figure 
3.4b: M e = 0.92, Figure 3.4c: M e = 0.80. 
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Figure 3.5. Full-field finite element results for the apparent stress exponent of 
the antisymmetric portion of the stresses at r/r, = 10' 6 and 10' 10 for a hardening 
exponent of n = 3 for plane strain. The results show the gradual transition from 
the HRR eigenvalue of -0.25 for mode II dominant loading to the eigenvalue S 2 
= -0.1275 for mode I dominant loading. These results confirm the asymptotic 
results in Figure 3.1. 
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Similar results have been obtained by Loghin (2001) for the case of plane stress. This work is currently 
being prepared for journal paper submission. For this case the full-field verification of the two different 
asymptotic crack solutions corresponding to Figure 3.2 are presented as Figures 3.6 and 3.7. 
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Figure 3.6. Comparison of the first and second terms of the mode II dominant 
asymptotic plane stress solution with the symmetric and antisymmetric parts of 
the ftill-field finite element solution at r/rj=10' 8 and r/r^lO" 4 for a hardening 
power of n = 3 and a mixed loading of M* = 0. 1 . 
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0.0 


- 0.1 
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Figure 3.7. Comparison of the first and second terms of the mode I dominant 
asymptotic plane stress solution with the symmetric and antisymmetric parts of 
the full-field finite element solution at r/ri=10' 8 for a hardening power of n = 3 
and a mixed loading of M e = 0.98. 
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This work is slow to be completed due to the rather confusing result of a singular power stronger 
than the HRR eigenvalue, for the antisymmetric part of the mode I dominant asymptotic response 
(see Figure 3.2). We believe this is explained by the fact that the asymptotic behavior is not 
actually valid in the mathematical limit as r approaches zero, in which case there is no violation of 
infinite energy. However, the result appears to be valid at extremely small values of r, well within 
the region where the HRR theory is no longer valid. To see this paradoxical behavior, in Figure 3.8 
results analogous to those of Figure 3.5 are presented for plane stress. The smallest value of (r/ri) 
where the theory can be applied is 10' 5 , which is on the order of the crack tip opening displacement. 



Figure 3.8. Plane stress case analogous to Figure 3.5. Apparent stress 
singularities from the full-field finite element analysis for a hardening 
exponent of n = 3 as a function of the loading parameter, Nf. The results at 
different values of ih\ are compared to the singularities given by asymptotic 
analysis in dashed lines. 


2. Higher order linear elastic terms 

A method to determine higher order terms from the numerical solution of a singular integral equation has 
been established. In addition to stress intensity factors and the T-stress, which are respectively the first and 
second terms, additional terms are determined. For the example of an edge crack or an internal crack in a 
half space, converged values of the first six coefficients have been determined, and for an edge crack in a 
finite width strip, the first four coefficients are given. Convergence of the results becomes more difficult for 
increasing numbers of terms. The results can be used to quantify the size of the K dominant zone and also 
the K-T dominant zone taking into account the entire region around the crack tip. The method can be used to 
determine cases where linear elastic fracture mechanics breaks down due to the loss of dominance of these 
zones. 


Williams (1952) was the first to express stresses and strains near the tip of a crack in terms of an asymptotic 
series for small distances from the crack tip. Using the polar coordinate system in Figure 1, the asymptotic 
form of the stresses and displacements for the symmetric, mode I case of loading can be expressed as 


a ij (r,0) = £kI(2r) n ^f i f(n,0) + XT:(2r)"f. In (n,0), 


(3.5) 


n=0 


n=0 
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(3-6) 


2^u i (r,e) = XkU2r) n+i g^(n,0) + |;T n , (2rr , g I , "(n,0),i = r,e,. 

n=0 n=0 

where the angular functions for both the symmetric and antisymmetric (mode II) cases are known. By 
applying considerable mathematical details and analysis, the coefficients can be very accurately determined 
using a singular integral equation approach. For example, for an edge crack in a half space the coefficients 
an infinite strip, the coefficients for tension and bending are given in Tables 3.1 and 3.2. 

Table 3.1 Asymptotic coefficients for an edge crack in a strip. 

Uniform tension. 


a/h 

K 

T 1 

A o 

k\b 

T,'b 

0.0 

a 0 ylb 

1.1215 

-0.5260 

0 o Vb 

0.2418 

°o 

-0.1925 

0.1 

1.1892 

-0.5502 

0.2460 

-0.19865 

0.2 

1.3673 

-0.5890 

0.21695 

-0.1994 

0.3 

1.1699 

-0.6103 

0.08426 

-0.1754 

0.4 

2.1114 

-0.57825 

-0.2873 

-0.09055 

0.5 

2.8246 

-0.4217 

-1.2293 

0.1617 

0.6 

4.0331 

0.03814 

-3.6968 

0.9509 

0.7 

6.3549 

1.3614 

-11.1150 

3.8845 

0.8 

11.9553 

6.0074 

-41.2088 

19.3852 

0.9 

34.6326 

35.7433 

-297.664 

213.725 

0.95 

99.1303 

166.299 

-1876.35 

1972.07 


Table 3.2 Asymptotic coefficients for an edge crack in a strip. 
Bending. 


a/h 

K 

T 1 

A o 

k[b 

T,'b 


cr,Vb 

al 

a,Vb 

a, 

0.0 

1.1215 

-0.5260 

0.2418 

-0.1925 

0.1 

1.0472 

-0.3779 

0.07967 

-0.08952 

0.2 

1.0553 

-0.2382 

-0.1112 

0.01899 

0.3 

1.1241 

-0.07917 

-0.3773 

0.14665 

0.4 

1.2606 

0.1208 

-0.7909 

0.3138 

0.5 

1.4972 

0.3975 

-1.51245 

0.5772 

0.6 

1.9140 

0.8339 

-2.9719 

1.1199 

0.7 

2.7252 

1.67525 

-6.6091 

2.6732 

0.8 

4.6764 

3.9269 

-19.3929 

9.5643 

0.9 

12.4621 

15.8057 

-116.047 

84.4986 

0.95 

34.3061 

63.2085 

-674.261 

713.191 
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For the case of an edge crack in a half-space, which is presented to show convergence capability, the first six 
coefficients are presented in Table 3.3. 

Table 3.3 Asymptotic coefficients for an edge crack in a half space. 


Term in Series 

klb" 

t y 

n 

a 0 Vb 


0 

1.1215222552 

-0.5259676011 

1 

0.241774599 

-0.1924915 

2 

0.027990 

0.0535? 


Section II: Measurables 

Publications this reporting period: 

1. Loghin, A. and Joseph, P.F., “Asymptotic Solutions for Mixed Mode Loading of Cracks and Wedges in 
Power Law Hardening Materials,” Engineering Fracture Mechanics. Vol. 68. pp. 1511-1534, (2001). 

2. Loghin, A. and Joseph, P.F., “Near mode I fracture in power law hardening materials,” International 
Journal of Fracture. Vol. 123, pp. 81-106, (2003). 

3. P.F. Joseph and M. Capitaneanu, “Determination of higher order coefficients from a singular integral 
equation approach,” in preparation. 

4. M. Capitaneanu, “Higher order terms in linear elastic fracture mechanics,” MS Thesis, Clemson 
University, May 2000. 

5. Loghin, A. and Joseph, P.F., “Mixed mode fracture in power law hardening material behavior for plane 
stress,” in preparation. 

6. A. Loghin, “The effects of arbitrary loading in nonlinear fracture mechanics,” Ph.D. Dissertation, 
Clemson University, August 2001. 

7. A. Aurora, P.F. Joseph, J.D. DesJardins, M. LaBerge, “The effect of lubricant composition on the fatigue 
properties of ultra high molecular weight polyethylene (UHMWPE) for total knee replacement (TKR),” 
submitted for publication, 2004. 

Grants and Financial Awards: 

Dr. Joseph has had one project initiate after the NASA EPSCoR project. This project is funded by the 

National Science Foundation through the NSF EPSCoR program. The total funding amount from the NSF is 

$499,000 over two years for five faculty. The research area is wear of polymers. Dr. Joseph’s research area 

in this NSF grant includes contact and fracture. 

Section III: Other Publications and Proposals: 

Presentations: 

1. P.F. Joseph, “Some complications in singular stress analysis,” University of Miami, Invited seminar, 
April, 2000. 

Section IV: Participants 

Graduate Students: 

Magda Capitaneanu (W,F) Earned MS degree in May 2000 

Full NASA funding, switched from Ph.D. program to MS. 

Adrian Loghin (W,M) Earned Ph.D. degree in August 2001 
Final year of Ph.D. program funded by NASA 

Amit Aurora (W,M) Earned MS degree in December 2003 
Funded by NASA for 6 months 
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Undergraduate Students: 

David Suttles (W,M) 

Faculty: 

Paul F. Joseph 

Section V: Collaborations 

NASA Centers: 

NASA Langley Research Center, Hampton, VA 
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(4) High Performance Titanium Alloys For Aerospace Applications 

H. J. Rack, Professor Materials Science and Engineering, Clemson University 

Dr. J. Qazi, Research Associate H.C. Bellum, M.S. Candidate 

I. Duta, M. S. Candidate S. Vignesh, M. S. Candidate 

Summary 

This task has examined a broad range of issues that limit the applicability of advanced titanium 
alloys for aerospace structures. Included are studies designed to define those factors that 
(a) control the deep hardenability and thus the heavy section capability of a + B and metastable B titanium 
alloys, (b)limit the high cycle fatigue behavior of metastable titanium alloys and (c) control the high 
temperature, high strain deformation behavior of a + B and metastable B titanium alloys. Finally the 
feasibility of designing a high modulus, high strength titanium alloy through coupled thermodynamic and 
solid mechanics modeling is being examined. 

Continued progress has been made in attaining the goals set forth in the original proposal. The study 
of phase transformations occurring during cooling of Ti-6A1-4V, IMI 550 and Ti-6-22-22 has been 
concluded. This have shown the importance of p phase stability in controlling type of martensite formed on 
rapid cooling, slower cooling resulting progressively in massive transformation, Widmanstatten, basket- 
weave and patchy a formation. Further the deep hardenability of metastable p titanium alloys has been 
shown to be limited by progressive nucleation of primary a at grain boundary triple points, grain boundary 
edges/surfaces and dislocation sub-boundaries - the critical cooling rates for each being dependent upon 
alloy composition. 

During this final reporting period studies of Ti-15V-3Cr-3Al-3Sn have shown that a 50% 
enhancement in high cycle fatigue performance can be achieved through controlled rolling/recrystallization 
and aging to form a bi-modal microstructure. Additionally it has been found that dynamic spheroidization in 
lamellae <x-p titanium alloys initially involves dynamic recovery/recrystallization of individual a lamellae 
followed by a/a grain boundary sliding and grain rotation, while flow stability in equiaxed a + B 
microstructures associated with dynamic recovery/recrystallization, unstable flow associated whith localized 
shear and strain assisted dissolution and coarsening of primary alpha phase. 

Finally coupled thermodynamic and solid mechanics modeling has been utilized to define a range of 
wrought Ti-Mo-Al/TiC alloys that have the potential of achieving both high modulus and high strength. 

Section I: Narrative 

Description of Research and Accomplishments 

The objective of this task is to establish an understanding of the factors that limit the development of 
advanced aerospace titanium alloys. These studies include examinations of 
phase transformations in a + B and metastable B titanium alloys, 
the high cycle fatigue behavior of metastable B titanium alloys, 
deformation processing of a + B and metastable B titanium alloys, and 
high modulus, thermally stable TiC reinforced titanium alloys. 

Phase Transformations 

This sub-task is presently examining phase transformations in a + B (IMI 550, Ti-6A1-4V) and 
metastable B (TIMET LCB and Beta-C) titanium alloys. 

a + P Alloys'. Previous year’s studies have shown that the transition from a hexagonal, a’ to an 
orthorombic, a”, martensite which occurs with increasing beta stabilizing content, and which limits the 
ability to attain high strength with adequate tensile ductility in solution treated and aged a + p titanium 
alloys, is related to the degree of short-range developed within the parent P phase, increasing degree of order 
promoting the formation of a” martensite. Qualitative agreement with this suggestion was achieved utilizing 
in-situ high temperature neutron diffraction at AECL, Chalk River, Canada. 
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Examination of the phase transformations taking place in EMI 550 during cooling from above or 
below the beta transus temperature as a function of cooling rate concluded. X-ray diffraction analysis of 
samples cooled from 1010°C (above the p-transus temperature) showed that the metastable P phase goes 
through the following transformation sequence with decreasing cooling rate, P~> a! (martensite) — > a! + a + 
P— » a + p. Microstructure analysis by optical and scanning electron microscopy of samples cooled at 
different rates showed that the samples also went through following microstructural changes; a! laths— » 
Widmanstatten a-» basket weave a-» basket weave + patchy a structure. The presence of a! martensite in 
the samples cooled from 1010°C at 450°C/sec detected by x-ray diffraction and optical microscopy was 
confirmed by transmission electron microscopy. The a! martensite laths in the samples cooled from 101 0°C 
at 450°C/sec has a highly dislocated structure, Figure 1, in agreement with the earlier work. X-ray diffraction 
analysis also revealed the presence of a and a! in the samples cooled at a rate in the range of 450 to 
100°C/sec, with the volume fraction of a 1 decreasing with slower cooling rates. Transmission electron 
microscopy of samples cooled at a rate of 100°C/sec confirmed the presence of a 1 and a, Figure 2; a being 
differentiated from a 1 by the presence of highly dislocated structure. Qualitative spot and line chemical 
analysis of these laths by STEM/EDX analysis showed partitioning of Mo and A1 in the a! and a laths 
respectively. Based on these results it is proposed that a' forms at M s temperature, further cooling resulting 
in auto-tempering of the martensite. 

X-ray diffraction analysis also revealed that on cooling from 900°C (below the beta transus 
temperature) IMI 550 goes through following phase transformations a pnmary + orthorhombic a! 1 martensite -> 
aprimaty + a 11 + P -» a pnmary + tXsecondaiy + P- Cooling from 900°C at a rate of 15°C/sec or faster resulted in the 
formation of orthorhombic a 11 martensite, with the volume fraction of a" decreasing with slower cooling 
rates. The microstructural changes were however too fine to be detected by optical or scanning electron 
microscopy, e.g., the a" martensite could not be detected by SEM analysis. Figure 3. The presence of a" 
martensite was confirmed by transmission electron microscopy, Figure 4. These fine a" needles have a 
highly twinned structure, in agreement with earlier work. 

The effect of cooling rate on the phase transformations in Ti-6A-4V on cooling from below the P- 
transus temperature has also been examined. Two temperatures, 940°C and 884°C, were selected to examine 
the role of primary a content and p phase composition on the phase transformation observed in this alloy. X- 
ray diffraction revealed that cooling from both these temperatures at 350°C/sec or faster resulted in the 
formation of a martensite. Additionally, XRD analysis showed that with slower cooling rates the P phase 
became stable and its volume fraction increased with decreasing cooling rate. TEM analysis of samples 
cooled at 350°C/sec from 940°C or 884°C confirmed the presence of a! martensite, having high dislocation 
density, stacking faults along with some twins. Additionally, TEM analysis of primary a in these samples 
confirmed its hexagonal structure its composition remaining unchanged at all cooling rates. Based on this it 
is suggested that for phase transformation purposes, primary a and the metastable p phases can be 
considered as two alloys, the former remains unchanged on cooling and later goes through a! to a+p 
transformation sequence with decreasing cooling rates. 

Metastable /3 Alloys: These studies during the final year of this program have focused on the transformations 
occurring during cooling of TIMET LCB and Beta-C. Transmission electron microscopy of TIMET LCB, as 
summarized in Figure 5, has shown that a phase nucleation initially occurs at grain boundary triple points, 
decreasing cooling rates resulting in progressive nucleation of a at grain boundary edges and surfaces, with 
homogeneous nucleation being observed at the slowest cooling rates examined. In contrast a nucleation 
during cooling in Beta-C occurs at substantially slower cooling rates with grain boundary, dislocation 
associated and homogeneous nucleation occurring progressively with decreasing cooling rate. 

High Cycle Fatigue Behavior of Metastable B Titanium Alloys: 

The previous results of this sub-task have shown that processing history can have a major impact on 
the high cycle fatigue performance of metastable p titanium alloys. Scanning electron microscopy has now 
shown that fatigue crack initiation is associated with the free surface for Ti-15V-3Cr-3Al-3Sn lot A sheet, 
Figure 6. In contrast, evidence of sub-surface initiation was found in lot B sheet, Figure 4. Fatigue crack 
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initiation in lot B also appeared to be associated with former a grain boundaries. Indeed, the transition in 
fracture appearance shown in Figure 7(insert) suggested that initiation in this sample occurred at a Tvpe 
I(coarse a)/Type II (fine a) interface. Examination of the fracture profile, also showed that fatigue crack 
propagation occurred predominantly through the Type I regions. 
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Figure 1. Transmission electron micrograph of IMI 550 solution treated at 1010°C followed by 
cooling at a rate of 450 3 C/sec showing (a) a laths, and (b) selected area diffraction pattern from 

these laths. 
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Figure 2. Scanning transmission electron micrograph showing a and a needles, carbon deposit on 
the laths show the location of EDX. 


Figure 3. Secondary scanning electron micrograph showing primary a particles in 1MI 550 sample 
cooled from 900°C at 150°C/sec. 
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Figure 4. Transmission electron micrograph showing a martensite needles in IMI 550 sample cooled from 

900°C at 150°C/sec. 
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Figure 5: Influence of cooling rate on a nucleation site in TIMET LCB. 

Nano-indentation measurements. Figure 8, of lot A specimens cycled at a max = 620 MPa, i.e., at the 
fatigue limit, demonstrated that failure at this stress level was associated with fatigue damage accumulation, 
that is cyclic hardening was observed in finite life fatigue samples, while cyclic softening was observed in 
samples where run-out occurred. Further measurements of lot B fatigue specimens, again cycled at the 
fatigue limit, a max = 827 MPa, indicated that the response in this material was dependent upon local 
microstructure. The cyclic response of Type I regions was similar to that of sheet A, while a continuous 
decrease in nano-hardness with increasing number of cycles was observed in Type II regions. These results 
indicate that enhancements in fatigue behavior within metastable (3 titanium alloys may be achieved through 
development of a bi-modal a microstructure. 

Deformation Modeling ; 

The effort during this past year within this task has focused on achieving an understanding of the 
high temperature, high strain deformation processes in lamella Ti-6A1-4V and equiaxed a+p Corona-X, the 
latter a higher toughness alloy developed in support of the HSCT program. Microstructural examination has 
shown that the flow behavior of lamellar Ti-6A1-4V under these conditions involves dislocation generation 
(from lamellae a/p interfaces), dynamic recovery and flow softening. Flow softening can be related to 
lamellae bending, kinking and dynamic spheroidization. The flow instability region defined by Dynamic 
Material Modeling at low temperatures and high strain rates has been related to the combined effect of shear 
localization and strain induced porosity at grain boundary aJ P interface. Finally electron back-scattered 
diffraction of individual a lamellae suggests that dynamic spheroidization involves dynamic 
recovery/recrystallization of a lamellae, followed by grain boundary sliding and rotation of these new a 
grains. 
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High temperature, high strain deformation of eqiuaxed microstractures has also been concluded. . 
These studies of Corona-X in the equiaxed a-p condition showed that stable flow, as defined by Dynamic 
Material Modeling techniques, can be associated with dynamic recovery and recrystallization, unstable flow 
being associated with localized shear(low' temperature and high strain rate), dynamic grain growth(high 
temperature and low' strain rate) and strain assisted dissolution and coarsening of the primary a phase(low 
strain rates at all temperatures examined). 



Figure 6: Scanning electron fractograph of Ti-15V-3Cr-3Al-3Sn lot A sheet in the solution-treated (ST) 
condition, aged condition ( 0 ^= 620 MPa, N= 65,689 cycles). 


I 
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Figure 7: Scanning electron fractograph of Ti-15V-3Cr-3Al-3Sn lot B sheet in the solution-treated (ST) 
condition, aged condition (a max = 827 MPa, N= 29,51 1 cycles). 
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Figure 8: Nano-indentation measurements performed on Ti-15V-3Cr-3Al-3Sn in the solution-treated (ST) 
condition, aged condition and fatigued specimens cycled at 0 ^= 620 MPa for lot A and a IQa x = 827 MPa for 

lot B sheet. 


Synthesis of TiC/Titanium Composites 

This sub-task examined the feasibility of designing a high modulus, high strength a + B titanium 
alloy utilizing particulate and dispersion strengthening techniques. Development of both Eshelby and finite 
element models of these three phase materials, containing a, B and TiC has been completed. These have 
shown that attainment of high moduli will depend upon control of carbide/alpha/beta phase size and 
distributions. The Ti-Mo-AFTiC system has been selected for this investigation. Thermodynamic 
simulations of this quaternary system using Thermo-Calc have shown that at elevated temperature 
approximately 1 wt pet C can be taken into solid solution, subsequent exposure at lower temperatures 
resulting in the precipitation of TiC dispersoids and reinforcements. Small buttons of the required 
compositions have been prepared and confirmed these predictions. 

Technology Transfer 

Phase Transformation in titanium alloys : Information on effect of cooling rate on phase 
transformations in IMI 550, Ti-6-22-22 and T-6A1-4V supplied to Boeing - Seattle to support F-22 program. 
Information on cooling rate effects on IMI 550 have been transmitted to GE Engine-Evandale, OH, that of 
TIMET LCB to TIMET. 

Fatigue Behavior of Metastable B Titanium Alloys : Information transmitted to TIMET. This 
information together with phase transformation studies on TIMET LCB being used to assist in 
implementation of TIMET LCB for automotive springs. 

Deformation Processing of Titanium : details of Ti-6A1-4V results provided to TIMET and ALLVAC 
to support internal processing activities. 

TiC Reinforced Titanium Alloys : Agreement signed with TIMET to provide support and scale-up of 
promising compositions as identified. Assistance has also been provided to Dynamet Technology, 
Burlington, MA to define optimized heat treatment schedule for TiC reinforced Ti-6A1~4V prepared from 
elemental powders. 

Honors and Awards 

National Materials Advisory Board, 2000 - 
Clemson University Board of Trustees Award for Faculty Excellence - 1997 

Humboldt Research Prize - 1997 
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Section II: Measurables 


Refereed Publications: 

1 . R. R. Boyer, H. J. Rack and V. Venkatesh, THE INFLUENCE OF THERMOMECHANICAL 
PROCESSING ON THE SMOOTH FATIGUE PROPERTIES OF Ti-15V-3Cr-3Al-3Sn, Mat'l Sci. Eng., 
Vol. A243, 7(1998). 

2. I. Philippart and H. J. Rack, HIGH TEMPERATURE DYNAMIC YIELDING IN BETA TITANIUM 
ALLOYS, Mat! Sci. Eng., Vol. A243, 196(1998). 

3. T. Ahmed and H. J. Rack, PHASE TRANSFORMATIONS DURING COOLING IN a + B TITANIUM 
ALLOYS, Mafl Sci. Eng., Vol.A243, 206(1998). 

4. I.Philippart and H. J. Rack, HIGH TEMPERATURE HIGH STRAIN DEFORMATION OF Ti-6.8Mo- 
4.5Fe-1.5Al, Mat’l Sci. Eng., Vol. A254, 253(1998). 

5. S. Azimzadeh and H. J. Rack, PHASE TRANSFORMATIONS IN Ti-6.8Mo-4.5Fe-l.5Al, Met. Trans., 
Vol. 29A, 2455(1998). 

6. M. Long, R. Crooks and H. J. Rack, HIGH-CYCLE FATIGUE BEHAVIOR OF SOLUTION 
TREATED METASTABLE B TITANIUM ALLOYS, Acta Met., Vol.47, 661(1999). 

7. S. Menon and H. J. Rack, FLOW STABILITY IN BINARY Al-Li ALLOYS, Mat’l Sci. Eng., Vol. 297, 
244(2001). 

8. S. Menon and H. J. Rack, HIGH TEMPERATURE, HIGH STRAIN FORGING BEHAVIOR OF 
BINARY Al-2Li, Mat’l Sci. & Tech., Vol. 17, 836(2001). 

9. K. Kharia and H. J. Rack, THE INFLUENCE OF SOLUTION TREATMENT TEMPERATURE ON 
THE MARTENSITIC PHASE TRANSFORMATIONS IN IMI 550(Ti-4Al-4Mo-2Sn-0.5Si), Met. 

Trans., Vol. 32A, 671(2001). 

10. A.Roussel, I. Duta and H. J. Rack, GRAIN GROWTH PHENEOMENA IN Ti-15V-3Al-3Cr-3Sn, 
Microstructural Modeling and Prediction during Thermomechical Processing. R . Srinivasan and S. L. 
Semiatin, eds., TMS, Warrendale, PA, 2001, pp. 157-165. 

1 1. A. Roussel, L. Reister and H. J. Rack, HIGH CYCLE FATIGUE RESPONSE OF SOLUTION 
TREATED AND AGED Ti-15V-3AI-3Cr-3Sn, FATIGUE 2002, Proc. 8 th Int. Fatigue Congress, 
Stockholm, Sweden, June 2002, A. F. Blom, ed., Vol. 3, pp. 1823-1832. 

12. M. Gheorghe, J.I. Qazi and H. J. Rack, ANISOTHERMAL ALPHA PHASE FORMATION IN Ti- 
6.8Mo-4.5Fe-1.5Al, Proc. 10 th World Titanium Conference, DGM, Hamburg, 2003, in press. 

13. J. Qazi, D. Hobson and H. J. Rack, PHASE TRANSFORMATIONS DURING COOLING OF Ti-6A1- 
2Mo-2Cr-2Sn-2Zr-0.12Si, Met. Trans., in preparation. 

14. J. Qazi and H. J. Rack, PHASE TRANSFORMATIONS DURING COOLING OF IMI 550, Met. Trans., 
in preparation. 

15. J. Qazi and H. J. Rack, PHASE TRANSFORMATIONS DURING COOLING OF SUB-TRANSUS 
SOLUTION TREATED Ti-6A1-4V, Met. Trans., in preparation. 

16. M. Gheorghe, J.I. Qazi and H. J. Rack, ANISOTHERMAL ALPHA PHASE FORMATION IN BETA - 
C TITANIUM, in preparation. 

17. T. Ahmed and H. J. Rack, PHASE TRANSFORMATIONS DURING COOLING OF Ti-5Al-2Sn-4Zr- 
2Cr-lFe(P-CEZ), Mat’ls Sci. & Tech., in preparation 

18. I. Duta and H. J. Rack, DYNAMIC SPHEROIDIZATION OF LAMELLAR TI-6A1-4V, in preparation.. 

19. S. Vignesh and H. J. Rack, High Temperature, High Strain Deformation of Equiaxed a+P CORONA -X, 
in preparation. 

Section III: Other Publications and Proposals 

Publications: 

1 . M. Long and H. J. Rack, TITANIUM ALLOYS IN TOTAL JOINT REPLACEMENT - A REVIEW, 
Biomaterials, Vol. 19,1621(1998). 

2. V. Venkatesh and H. J. Rack, ELEVATED TEMPERATURE STRAIN HARDENING OF INCONEL 
690, Mech. of Mat’ls, Vol. 30, 69(1998). 
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3 . V. Venkatesh and H. J. Rack, INFLUENCE OF MICROSTRUCTURAL INSTABILITIES ON 
ELEVATED CREEP DEFORMATION OF INCONEL 690, MatTs Sci. & Tech., Vol.15, 408(1999). 

4. V. Venkatesh and H. J. Rack, A NEURAL NETWORK APPROACH TO ELEVATED 
TEMPERATURE CREEP-FATIGUE LIFE PREDICTION, Int. Jn Fatigue, Vol. 21, 225(1999). 

5. A. C. Geiculescu, H. G. Spencer, H. J. Rack and B. Sullivan, AQUEOUS SOL-GEL MODIFICATION 
OF THE FIBER-MATRIX INTERFACE IN GRAPHITE FIBER ALUMINUM MATRIX 
COMPOSITES, Adv. Mat'ls & Man., Vol. 14, 489(1999). 

6. X. Tang, T. Ahmed and H. J. Rack, MARTENSITIC PHASE TRANSFORMATION IN Ti-Nb-Ta and 
Ti-Nb-Ta-Zr ALLOYS, Jn. Mat'ls Sci., Vol. 35, 1805(2000). 

7. J. Fitzgerald, P. D. Rack, A. C. Geiculescu and H. J. Rack, DEPOSITION OF CERAMIC MATERIALS 
USING POWDER AND PRECURSOR VEHICLES VIA DIRECT WRITE PROCESSING, Materials 
Development for Direct Write Technologies. MRS Proc. Vol. 624, 2000, pp. 29-34. 

8. A. C. Geiculescu and H. J. Rack, ATOMIC-SCALE STRUCTURE OF WATER-BASED XEROGELS 
BY X-RAY DIFFRACTION, Jn. Sol-Gel Sci. & Tech., Vol. 30, 13(2001). 

9. P. D. Rack, A. C. Geiculescu and H. J. Rack, GROWTH OF REACTIVE SPUTTER DEPOSITED 
TUNGSTEN-CARBON THIN FILMS, Jn. Vac. Sci. Tech., Vol. 19, 62(2001). 

10. M. Long and H. J. Rack, FRICTION AND SURFACE BEHAVIOR OF SELECTED TITANIUM 
ALLOYS DURING RECIPROCATING-SLIDING MOTION, Wear, Vol.249, 158(2001). 

1 1. A. C. Geiculescu and H. J. Rack, ATOMIC SCALE STRUCTURE OF ALCOHOL-BASED ZIRCONIA 
XEROGELS BY X-RAY DIFFRACTION, Jn. Non-Ciystalline Solids, Vol. 289 53(2001). 

12. H. J. Rack, I. Gheorghe, K. Kharia and A.C. Geiculescu, IN-SITU FABRICATION OF FIBER 
REINFORCED METAL MATRIX COMPOSITES, Affordable Metal Matrix Composites for 
High Performance Applications. A. Pendey, K. Kendig and T. Watson, eds., TMS, Warrendale, 
PA, 2001, pp. 21 1-221. 

13. I.Gheorghe and H. J. Rack, INFLUENCE OF Ti0 2 TO A1 RATIO ON THE REACTION IN Ti0 2 /Al 
COMPOSITES, MatTs Sci and Tech., Vol. 18, 1079(2002). 

14 . I.Gheorghe and H. J. Rack, REACTIVE INFILTRATION OF 25 V/O Ti0 2 /Al COMPOSITES, Met. 
Trans., Vol. 33 A, 2155(2002). 

15. A. C. Geiculescu and H. J. Rack, SMALL ANGLE X-RAY SCATTERING STUDIES OF 
POLYMERIC ZIRCONIA AQUEOUS XEROGELS, Jn. Non-Crystalline Solids, Vol. 306, 30(2002). 

16. G.G. Yapici, I. Karaman,Z. P. Luo and H. Rack, MICROSTRUCTURE AND MECHANICAL 
PROPERTIES OF SEVERELY DEFORMED POWDER PROCESSED Ti-6A1-4V USING 
EQUAL CHANNEL ANGULAR EXTRUSION, Scripta Materialia, Vol. 49,1021(2003). 

20. M. Long, J.I. Qazi and H.J. Rack, RECIPROCATING SLIDING WEAR RESISTANCE OF |321SRx 
TITANIUM, Proc. 10 th World Titanium Conference, DGM, Hamburg, 2003, in press. 

21. J.I. Qazi, V. Tsakiris, B. Marquardt and H. J. Rack, THE EFFECT OF DUPLEX AGING ON THE 
TENSILE BEHAVIOR OF Ti-35Nb-7Zr-5Ta-(0.06-0.7)0 ALLOYS, Proc. 10 th World Titanium 
Conference, DGM, Hamburg, 2003, in press. 

22. J. I. Qazi and H. J. Rack, EFFECTS OF THERMAL TREATMENT ON THE MECHANICAL 
PROPERTIES OF BIOMEDICAL TITANIUM ALLOYS, MatTs and Processes for Medical 
Devices I. ASM Int., in press. 

23. R. Z. Valiev, V. V. Stolyarov, H. J. Rack and T. C. Lowe, SPD-PROCESSED ULTRA-FINE 
GRAINED TI MATERIALS FOR MEDICAL APPLICATIONS . MatTs and Processes for 
Medical Devices I . ASM Int., in press. 

24. M. Kocan, T. Ludian, M. Ishii, H, J, Rack and L. Wagner, OPTIMIZATION OF MICROSTRUCTURE 
OF TIMET AL LCB FOR APPLICATION AS SUSPENSION SPRINGS, Light Materials for 
Transportation SvstemsCLIMATL ed. O. Es-Said, 2003, in press. 

25. E. Fu, H. C. Bellum, J. Qazi, H. J. Rack and V. Stolyaru, RECIPROCATING-SLIDING WEAR OF 
ULTRA-FINE GRAINED Ti-6A1-4V, Proc. 3 rd Int. Sym. On Ultrafine Grained Materials, TMS, in press. 

26. M. Long and H. J. Rack, RECIPROCATING-SLIDING WEAR AND SUBSURFACE BEHAVIOR OF 
TITANIUM ALLOYS, Met. Trans., in press. 
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27, A. C. Geiculescu, K. Siegamura and H. J. Rack, REACTIVE INFILTRATION OF 25 V/O Nb 3 0 5 /Al 
COMPOSITES, in preparation. 

Presentations: 

1 . H. J. Rack, PHASE TRANSFORMATIONS IN METASTABLE BETA TITANIUM ALLOYS, 

Aeromat ’97, Williamsburg, VA, May, 1997. 

2. H. J. Rack, HIGH TEMPERATURE DEFORMATION PROCESSING OF ct+P TITANIUM ALLOYS, 
Aeromat ’97, Williamsburg, VA, May 1997. 

3. R. Boyer, H. J. Rack and V. Venkatesh, THE INFLUENCE OF THERMOMECHICAL PROCESSING 
ON THE SMOOTH FATIGUE PROPERTIES OF Ti-15V-3Cr-3Al-3Sn, International Conference on 
Thermomechanical Processing of Steels and Other Materials, Univ. of Wollengong, Australia, July, 1997. 

4. T. Ahmed and H. J. Rack, PHASE TRANSFORMATIONS DURING COOLING IN a+p TITANIUM 
ALLOYS, International Conference on Thermomechanical Processing of Steels and Other Materials, Univ. 
of Wollengong, Australia, July, 1997. 

5. H. J. Rack and I. Phillipart, HIGH TEMPERATURE DYNAMIC YIELDING IN BETA TITANIUM 
ALLOYS, International Conference on Thermomechanical Processing of Steels and Other Materials, Univ. 
of Wollengong, Australia, July, 1997. 

6. J. Brasell, P. Martin and H. J. Rack, PHASE TRANSFORMATIONS IN CORONA-X, Fall Mtg 
ASM/TMS, Indianapolis, Indiana, Sept. 1997. 

7. X. Tang and H. J. Rack, KINETICS OF ORDERING IN Ti-6Al-2Mo-2Cr-2Zr, Fall Mtg ASM/TMS, 
Indianapolis, Indiana, Sept. 1997. 

8. H. J. Rack, PRECIPITATION IN METASTABLE BETA TITANIUM ALLOYS, Tech. Univ. Hamburg- 
Harburg, June 1998. 

9. K. Kharia and H. J. Rack, PHASE TRANSFORMATION IN IML 550, TMS Annual Meeting, Nashville, 
Tenn., March, 2000. 

10. H. J. Rack and C. Robinson, HIGH TEMPERATURE, HIGH STRAIN DEFORMATION OF BETA-CEZ 
TITANIUM, TMS Annual Meeting, Nashville, Tenn., March 2000. 

11.1. Gheorghe and H. J. Rack, THERMODYNAMICS OF IN-SITU REACTIONS BETWEEN Ti0 2 AND 
PURE AL, Annual Mtg, TMS, Nashville, Tenn., March, 2000. 

12. A. C. Geiculescu and H. J. Rack, SHORT RANGE STRUCTURE OF ZIRCONLA XEROGELS, Am. 
Cer. Soc., Annual Mtg., St. Louis, Mo., April 2000. 

13. J. Fitzgerald, P. D. Rack, A. C. Geiculescu and H. J. Rack, DEPOSITION OF CERAMIC MATERIALS 
USING POWDER AND PRECURSOR VEHICLES VIA DIRECT WRITE PROCESSING, MRS Mtg., 
San Francisco, CA, April, 2000. 

14. H. J. Rack, ADVANCED METALLIC MATERIALS FOR BIOLOGICAL APPLICATIONS, INSA- 
Rennes, France, June 7, 2000 and TUMunich, Germany July 20, 2000. 

15. H. J. Rack, HIGH TEMPERATURE, HIGH STRAIN DEFORMATION OF NEAR-BETA 
TITANIUM ALLOYS, BTU-Cottbus, Germany, June 29, 2000. 

1 6 . H. J. Rack, MARTENSITIC PHASE TRANSFORMATIONS IN ALPHA-BETA TITANIUM ALLOYS, 
BTU-Cottbus, Germany, July 3, 2000. 

1 7 . H. J. Rack, IN-SITU FORMATION OF HIGH TEMPERATURE COMPOSITES, BTU-Cottbus, 
Germany, July 5, 2000. 

1 8 . H. J. Rack, IN-SITU FORMATION OF LIGHT- WEIGHT FIBER REINFORCED COMPOSITES, ARL, 
Belvedere, Md, August 21, 2000. 

19. H. J. Rack, J. Wendt, M. Hilpert and L. Wagner, DIRECTIONAL MECHANICAL PERFORMANCE OF 
WROUGHT AZ80 MAGNESIUM, Annual Mtg, TMS, New Orleans, Feb. 17, 2001 . 

20. H. J. Rack, TITANIUM IN THE 21 st CENTURY, BTU-Cottbus, Germany, October 26, 2001. H. J. Rack, 
ANISOTHERMAL ALPHA PHASE FORMATION IN Ti-6.8Mo-4.5Fe-l.5Al, 10 th World Titanium 
Conference, DGM, Hamburg, July 11, 2003. 

21. H.J. Rack, RECIPROCATING SLIDING WEAR RESISTANCE OF p21SRx TITANIUM, 10 th World 
Titanium Conference, DGM, Hamburg, July 12, 2003. 
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22. H. J. Rack, THE EFFECT OF DUPLEX AGING ON THE TENSILE BEHAVIOR OF Ti-35Nb-7Zr- 
5Ta-(0.06-0.7)0 ALLOYS, 10* World Titanium Conference, DGM, Hamburg, July 12, 2003. 

23. H. J. Rack, EFFECTS OF THERMAL TREATMENT ON THE MECHICAL PROPERTIES OF 
BIOMEDICAL TITANIUM ALLOYS, Mat’ls and Processes for Medical Devices, Anaheim, Cal., Sept. 
8, 2003. 

24. H. J. Rack, SPD-PROCESSED ULTRA-FINE GRAINED TI MATERIALS FOR MEDICAL 
APPLICATIONS, Mat’ls and Processes for Medical Devices, Anaheim, CA, Sept. 8, 2003. 

Proposals: 

1 . Phase Transformations in TiC Reinforced P/M Ti-6A1-4V, DOE, March 2002 

2. Interrelationship between Microstructure and High Cycle Fatigue Performance of TIMET LCB, TIMET, 
July 2002(Funded). 

3. Novel Titanium Materials for Orthopaedic Lumbar Implants, NTH, December, 2002(Funded). 

Section IV: Participants 

Current Participants 

Graduate Students: 

M. S. Candidates: 

Ion Duta (M,W) 

S. Vignesh(M, Asian) 

H. C. Bellam(F, Asian) 

Research Associates: 

Dr. J. Qazi(M, Asian) 

Faculty: 

Henry J. Rack (W, M) 

Former Participants 

Undergraduate Students: 

Rebecca. Hubbard (F, W) 

Graduate Students: 

M. S. Graduates: 

John Brasell (M, W)[Nucor Steel, Darlington, SC] 

Alex Rousell (M, W)[ALLVAC, Monroe, NC.] 

Dale Hobson (M, W)[Oxford Instruments, Medford, Mass] 

Krishna Kharia (M, Asian)[Tata Industries, India] 

Mariana Gheorghe (F, W)(PhD Candidate Georgia Tech.) 

Ph.D. Candidates: 

Santosh. Menon (M, Asian)[LSI Logic, Portland, OR] 

Research Associates: 

Dr. T. Ahmed (M, Asian), SABIC, Riyadh, Saudi Arabia 
Dr. X. Tang (F, Asian)[Carpenter Technology, Reading, Pa] 

Section V: Collaborations 

NASA Centers: 

NASA Langley Research Center. Hampton. VA: 

Mr. William Brewer - HSCT Program 
Mr. Dennis Discus - HSCT Program 

Ms. Terry! Wallace, serve on PhD committee at Univ. of Virgina 

Dr. Roy Crooks, joint research on grain boundary misorientation effects on high 

cycle fatigue behavior of metastable beta titanium alloys 
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Academic: 


Tech. University of Munich, Germany 

Dr. Jean Gregory, high cycle fatigue performance of metastable P 
titanium alloys 

BTU- Cottbus/TUClausthal , Germany 

Dr. Lothar Wagner, fatigue and fracture of titanium alloys 
INSA-Rennes, France 

Dr. Jean Debuigne, phase transformations in metastable P titanium 
alloys 

Business and Industry: 

Titanium Manufacturers: 

TIMET, Las Vegas, NV : 

Dr. Steve Fox, Design of TiC/Ti Composites, fatigue performance of 
metastable beta titanium alloys. 

Dr. Vasisht Venkatesh, former graduate student 
RMI, Niles, Ohio: 

Mr. James Ferraro, processing of EMI 550, Ti-6-22-22, Beta-C 
Teledyne ALLVAC, Monroe, NC: 

Mr. Alex Roussel (former student), Dr. Richard Kennedy, Mr. J. R. Wood, 
Mr. J. R. Wood, Mr. Michael Volas, processing of titanium alloys 
Airframe Manufacturers: 

Boeing Commercial Airplane, Seattle, WA: 

Dr. Jim Cotton, Mr. Rodney Boyer, thermomechanical processing of 
metastable beta titanium alloys 
Boeing, St. Louis, MO: 

Mr. Richard Lederich, phase transformations in Ti-6-22-22 
Aircraft Engine Manufacturers: 

General Electric, Evandale, Ohio: 

Dr. Andy Woodfield, processing of a-B titanium alloys. 

Tom Broderick, processing of IMI 550 
Other Federal/Governmental Agencies: 

Oak Ridge National Laboratories. High Flux Reactor Site: 

Dr. Stephan Spooner, Neutron Diffraction Analysis of Ordering in IMI 550 
AECL. Chalk River. Canada 

Dr. Kelly Conlon, High Temperature Neutron Diffraction 
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NASA EPSCoR Final Report 


Name of PI - Nikunja K. Swain, Ph.D., P.E. 
Cluster - Aeronautics 


Section I: Narrative 


I worked on this project during summer 2003 for only six weeks. My research focused on 
digital signal processing using Joint Time and Frequency Analysis technique. 

Need for JTFA 

A given signal can be represented either as a function of time, which shows how the 
signal amplitude changes over time or as a function of frequency, which tells us how 
frequently the amplitude changes. The bridge between time and frequency 
representations is the Fourier transform, which is one of the most important signal 
analysis methods. Using Fourier transformation, a signal can be easily decomposed as a 
weighted sum of sinusoid functions. This enables one to process either the signal time 
waveform or its corresponding set of sinusoid functions, depending on which form is 
more convenient. Fourier transform also used to compute power spectrum for a signal. 
The power spectrum usually has a simpler pattern than the time waveform, and as result it 
often serves as the fingerprint of the analyzed signal. But the Fourier transform possesses 
certain disadvantages that prevent its use to represent many signals encountered in real 
world effectively. Examples of some of these signals are seismic and ECG 
(electrocardiogram), speech and noise corrupted signals. To properly analyze this type of 
signals both time and frequency analysis are generally performed. This is called joint 
time-frequency analysis (JTFA). The main objective of this analysis is to determine the 
energy density of a signal, which indicates how much energy is contained within a certain 
time interval in a given frequency band. This information is then used to analyze the 
signal properties. 


I continued designing and developing Virtual Instrument instructional modules for digital 
signal processing and digital filters. I also spend time on design and development of 
remote control and data acquisition modules using Lab VIEW and Visual Basic (Object 
Oriented Programming). 

I have prepared a manuscript entitled “Joint Time Frequency Analysis Using Virtual 
Instruments” and plan to present it in local/regional conference. 

The project budget provided funding only for my summer salary and there was no set 
aside funding in my part of the proposal for involving students. Therefore, no student was 
involved in my part of the project. 
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Section II: Measurables 


Publications - 

1. N. K. Swain, J. A. Anderson, R. Korrapati, M. Swain, Ajit Singh, “IPAR - A Visual Basic 

Program to study IP Addressing and Routing” - Proceedings of the IEEE SECoN 2003, 
Jamaica, West Indies, April 2003. 

2. N. K. Swain, Ajit Singh, J. A. Anderson, Mswain, et. all, “Remote Data Acquisition, 

Control, and Analysis using Lab VIEW Front Panel and Real Time Engine” - Proceedings of 
the IEEE SECoN 2003, Jamaica, West Indies, April 2003. 

3. N.K. Swain, Mrutyunjaya Swain, James A. Anderson, “Integrating Virtual Instruments into 
EET Curriculum - 2004 ASEE Annual Conference Proceedings, Salt Lake City, Utah, 

June 2004. 

Patents - None 

Grants and Financial Awards - 

“Computer Based Virtual Engineering Laboratory - CBVEL” - DOD Infrastructure grant, 1999. 

Section III - Other Publications and Proposals 

None 

Section IV - Participants 
None 

Section V - Collaborations 


None 
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Section i: Executive Summary of 
the 5 th and Final Year of NASA 
Sponsored EPSCoR Research in 
South Carolina 

This report summarizes both the Wetlands 
Research Cluster’s Year 5 activities and our 
research accomplishments over the full five 
years of our grant. This grant has been a 
tremendous asset to the colleges involved and 
a tremendous boost to environmental research 
within South Carolina. 

Quite notably, our remote sensing cluster has 
matured through our NASA EPSCoR 
collaborative research and experiences. True 
to the program’s intent, several of our team 
members have won additional competitive 
research awards totaling more than $4.5M 
from 6 non-NASA agencies an have additional 
proposals pending for an additional $5M. All 
members of our team were active in outreach: 
more than 30 presentations were delivered to 
local, state, regional and national audiences, 
and we published more than 50 papers and 
abstracts. We are very proud of the thirty-three 
M.S. and Ph.D. students we supported with 
this funding. Many of these students are now 
Pi’s and are now writing their own research 
grants. We also supported 9 undergraduate 
students who are now in graduate school or 
recently graduated from environmental 
science programs. Two student internship 
opportunities were also initiated; one 
corporate and one with SC Department of 
Natural Resources. 

In addition to gathering and analyzing a variety 
of remote sensor data to evaluate the health of 
South Carolina’s lowcountry wetlands, our 
team also looked at the effects of population 
growth on the environment. A full report of the 
culmination of 5 year’s worth of modeling 
population growth is presented in Appendix 1 . 

Year 5 & Final Research Results: 
Incorporation of BRDF Characteristics of 
Spartina alterniflora to Model Estuarine- 
wide Biomass Characteristics and 
Perform a Change Detection Analysis. 

Building upon our advancements in the use of 
innovative remote sensing and digital image 
processing techniques over the past four 


years, we performed a change detection 
analysis of biomass characteristics of Spartina 
alterniflora within the Murrells Inlet estuary. 
During the summer of 2003, we acquired high 
resolution multi-spectral data of the Murrells 
Inlet estuary using the same protocol for the 
acquisition of the remote sensing and in situ 
calibration data as developed during Year 1 of 
this project (Jensen et al., 1998). 


Activity 1: Murrells Inlet 1997-2001 
Change Detection Analysis 

We quantified and investigated the changes in 
upland and vegetative cover in Murrells Inlet, 
SC over the five-year grant period. A change 
detection analysis will be conducted using 
ATLAS 3 x 3 m multispectral data acquired on 
October 7, 1997 and IKONOS 4 x 4 m multi- 
spectral data acquired in the fall of 2001. 
Potential anthropogenic-induced trends (i.e., 
rates of increases in impervious surface) were 
identified that demonstrate the correlations 
between estuarine vegetation stress and 
upland land cover changes. Research 
previously conducted of BRDF signatures of 
Spartina alterniflora were incorporated into the 
ATLAS and IKONOS data for improving image 
classification. Radiometric adjustments that 
account for solar and sensor zenith angles 
were used for deriving accurate surface 
albedos. 

In addition, in-situ estuarine biophysical data 
(i.e., biomass, LAI, height, etc.) were acquired 
simultaneously with the IKONOS data 
collection for modeling changes in the 
biophysical parameters of the estuary between 
1997 and 2001. 


Activity 2: Evaluating ADAR or Hymap 
data for mapping intertidal oyster reefs. 

The eastern oyster, Crassostrea virginica, 
serves prominent roles in the ecology of 
estuarine and coastal ecosystems (Bushek, et 
al., in press). This oyster, which lives in the 
intertidal portion of South Carolina’s estuaries, 
is a reef-forming, filter feeding bivalve. Oyster 
reefs perform a variety of ecological roles that 
contribute significantly to the structure and 
function of the coastal ecosystems in which 
they exist. Oyster reefs form a unique habitat 
that supports and abundance of taxonomically 
diverse species. They stabilize creek banks 
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and serve as both active and passive filters as 
waters pass over them. In South Carolina, the 
eastern oyster is both an important for 
commercial fishery and recreational fishery 
activities. The natural population of oysters 
has declined dramatically due to overfishing 
and habitat loss/degradation (Bushek 2001). 

In South Carolina, the SC Department of 
Health and Environmental Control is 
responsible for regulating the opening and 
closing of shellfish reefs for commercial and 
recreational harvesting. The SC Department of 
Natural Resources (SCDNR) is responsible for 
managing state shellfish ground for 
commercial and recreational harvesting. 
Important to the efforts of both agencies is the 
ability to accurately map existing shellfish 
resources in a timely manner. Currently, the 
SCDNR maps shellfish reefs using a 
combination of Global Positioning System 
(GPS) surveys and field measurements and 
incorporates the information in a Geographic 
Information System (GIS). While much 
improved over historical techniques which 
involved extensive field surveys and laborious 
manual cartographic techniques are still 
manpower intensive and limit the ability of the 
agency to maintain and update oyster maps 
for the entire coast of South Carolina. 

Therefore, the purpose of this study is to 
assess the utility of supervised classification of 
high spectral resolution imagery for mapping 
oyster beds in a South Carolina estuary. The 
proposed study site is the intertidal area of the 
North Inlet-Winyah Bay National Estuarine 
Research Reserve located north of 
Georgetown, SC. The imagery used in this 
study is a rectified and registered image 
acquired by Positive Systems, Inc. at low tide 
on November 4, 1999 and HYMAP 

hyperspectral imagery obtained on October 
22, 2000. The ADAR 5500 digital camera was 
used to acquire 4 bands (band 1 : 450-51 5 nm; 
band 2: 525-605 nm; band 3:630-690 nm; and 
band 4: 750-860 nm) at 0.7 x 0.7 m spatial 
resolution. The imagery was acquired in 
support of this project and an associated 
Environmental Protection Agency project. 

Building upon a student project (Schmidt 
2000) supported by this NASA grant, a GER 
1500 spectroradiometer was used to acquire 
multiple spectral profile samples of oyster 
beds throughout the study area. Profiles of live 


and empty oyster shells, pure mud, mixed with 
dead oyster shells will be acquired. Profile 
information, logistic regression techniques, 
and image processing techniques, were used 
to develop a supervised classification of 
shellfish beds within the intertidal zone of the 
North Inlet study area. We also attempted to 
classify the health (i.e., living vs. dead) of the 
oyster beds using these techniques. The 
resultant map was ground-referenced to 
determine the accuracy of this use of remotely 
sensed imagery for mapping shellfish 
resources. 

As part of this project, we established two 
student internship opportunities. One is with 
Geometries, an engineering consulting firm 
specializing in coastal environmental issues, 
and the second with the Marine Resources 
Division of the SCDNR. 

References 

Bushek, D. 2001 . 

http://marinescience.sc.edu/facultv/bushek.sht 
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http://www.cla.sc.edu/aeoq/rslab/751/f00/schm 

idt/ 

Bushek, D., J. Kreesee, B. Jones, D. White, 
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health management: a system level 

perspective for Perkinsus marinus. Journal of 
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Activity 3: Evaluating the use of 
HYMAP imagery for mapping algal 
community composition and wetland 
biomass. 

Qualification and quantification of algal 
community composition in estuarine systems 
can yield important information about overall 
water quality, eutrophication, and toxic and 
non-toxic algal blooms. Common sampling 
techniques have relied upon point sampling 
from land and watercraft sampling stations, 
and data analyses have depended upon highly 
trained individuals using microscopy. Although 
these techniques are effective, they are limited 
in spatial coverage and the time and effort 
required to analyze a large number of samples 
is prohibitive. 
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Recent advances in the use of high 
performance liquid chromatography (HPLC) 
have aided researchers to qualify algal 
community composition (Wright et al., 1996). 
Further research in the use of factor analysis 
to analyze HPLC results has led to a robust 
method to quantify each component of the 
algal community can be described by 
accessory pigment(s) or pigment ratios in 
relation to total chlorophyll a (Mackey et al, 
1 996). Thus, if the pigment characteristics of 
the local algal communities are adequately 
described from lab cultures, wild algal 

populations from estuaries along the South 

Carolina coast (White and Lewitus 

unpublished data). 

The use of remotely sensed data from aircraft 
and spacecraft to quantify surface water algal 
biomass by monitoring chlorophyll a has 
become common, allowing for large spatial 
coverage repeatable over time. However, the 
available remote sensing platforms were 

limited to chlorophyll a detection. The advent 
of hyperspectral data may hold promise in the 
use of remotely sensed data may hold 
promise in the use of remotely sensed data to 
derive an overall pigment composition from a 
water body. The Airborne Visible InfraRed 
Imaging Spectrometer (AVIRIS) was the first 
of these platforms to be used to classify and 
map a mixed phytoplankton bloom 
(Richardson and Kruse). 

The purpose of this study is to assess the 
HYMAP sensor (Hyperspectral Mapping) in 
describing surface water algal community 
composition in a well-mixed coastal estuarine 
system. The HYMAP sensor was deployed 
over North Inlet, South Carolina in October of 
2000. This sensor has a range of 100 - 200 
bands and a bandwidth of 10-20nm. In 
conjunction with this over flight, water samples 
were collected and a suite of environmental 
variables were analyzed including inorganic 
and organic nitrogen and phosphorus, algal 
pigments ( chi a and accessory pigments), total 
sediment load and total organic carbon. A 
GER 1500 spectro-radiometer was used to 
acquire multiple spectral profile samples of the 
water surface from sampling stations. The in- 
situ spectroradiometer profile information and 
specialized hyperspectral digital image 
processing techniques were used to classify 
the imagery. The techniques included 
matched filtering, spectral feature fitting, and 


pixel purity index analysis. The hyperspectral 
data will also be analyzed and compared to 
HPLC and microscopy results to determine the 
accuracy of the remotely sensed data. If 
possible, the hyperspectral data were 
evaluated to determine the biomass of the 
smooth cordgrass ( Spartina altemiflora) 

surrounding the estuarine water under 
investigation. Results are discussed more 
fully in Section II. 


HYMAP Sensor Specifics 
HyMAP sensor is a 126-band commercial 
hyperspectral scanner that is capable of 
acquiring data at spatial resolution between 3- 
10m. The spectral resolution is listed below: 
VIS 0.45 - 0.89 um (15 nm sampling interval) 
NIR 0.89 - 1 .35 um (15nm sampling interval) 
SWIR1 1.40-1.80 um (13 nm sampling 
interval) SWIR2 1.95 - 2.48 um (17 nm 
sampling interval) 

Additional information: 

http://www.analvticalimaaina.com/data aco/hv 
map 1999/hvm specs.htm 

HvMAP Data report for 10-22-2000 North 
Inlet. SC 

Geo: average alt: 2098.685895 

Geo: Start line lat, Ion, alt: 33.381052, - 

79.170670, 2087.392954 

Geo: End line lat, Ion, alt: 33.297326, - 

79.214087, 2102.480411 

Geo: year, month, day: 2000, 10, 22 

Geo: hour, minute, sec: 16, 17, 40 

Geo: average scene elevation: 2.233238 

Geo: UTM zone: 17, Datum: NAD-27 

Geo: original pixel size: 4.500000 

Geo: output pixel size: 4.00000 

Geo: UTM upper left X coord (1,1) 

665081.445414 

Geo: UTM upper left Y coord (1,1) 
3695302.511934 
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Figure 1. Murrells Inlet land cover for 1997 and 2001. 


Section II: Measurables for 
Wetlands Remote Sensing 
Activities 

Research efforts towards quantifying upland 
and vegetation cover in Murrells Inlet, SC 
succeeded in data acquisition and remote 
sensing data analysis. Changes in sensor 
availability and specific requests made by 
local resource management organizations, 
however, altered the choice of airborne sensor 
and specific tasks, though the original 
research requirements remained. 


In support of a request from the Murrells Inlet 
2005 Task Force and the NOAA-supported 
Murrells Inlet Special Area ManagementPlan 
(MI-SAMP), an effort was undertaken to 
provide a land cover change detection 
analysis using imagery collected as part of this 
project and an additional NOAA-funded 
research effort. Using CAMS imagery acquired 
in August 1997 in support of this effort, and 
IKONOS imagery acquired in October 2001 in 
support of the NOAA-funded project, a change 
detection analysis was performed. The results 
of the change detection analysis suggest the 
most significant change within the Murrells 
Inlet area was an increase in impervious 
surfaces (Figures 1 and 2). 
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Area in Acres 



Shrub Wettand Impervious Sand Water Managed Grass Forest Barren 

Land Cover 

Figure 2. Results of 1997 to 2001 land cover change detection analysis. 


As part of an ongoing coast-wide initiative to 
improve the ability to map and assess 
intertidal shellfish resources, the South 
Carolina Department of Natural Resources 
requested that we evaluate imagery acquired 
as part of this effort to determine utility for 
mapping the spatial distribution and 
characterizing the structure of mapped 
shellfish beds. 

The selection of study sites was designed with 
the idea of gathering spectral signatures of 
shellfish in various environmental conditions. 
Conditions such as: abundance of shellfish, 
positions, (vertical and horizontal), and 
environmental condition, (wet/dry, various 
amount of mud present, algae, and detritus). 
Site characterization of oyster reefs utilized 
definitions from the Shellfish Management 
section of the South Carolina Department of 
Natural Resource’s Intertidal Oyster Survey 
Field Data Sheet codes. These describe the 
reef strata with respect to bushels of live 
oysters per acre, presence or absence of 
vertical clusters, proportion of live oysters to 
shells and amount of mud present. Each of the 
1 0 listed strata or classifications are identified 
by a letter code. Field collection using the 
hand-held GER spectroradiometer to gather 
spectral signatures of shellfish, vegetation, 


mud, and water took place over a 12-month 
period and constitute the core of the spectral 
library that will be used to derive spectral 
endmembers from remotely sensed imagery. 
Figure 3 showing the sample site, BOB4. 
Twenty-two sample points in total were 
established on BOB4 representing the various 
shellfish strata’s present. Figure 3 shows five 
of those sample points and the next four 
figures illustrate the spectral signatures of 
mud, shellfish, mud and shellfish, and a 
combined chart of all spectral signatures 
including Spartina altemiflora as percent 
reflectances. 

By comparing the percent reflectances otmud, 
shellfish and vegetation we see there is a 
difference between mud and shellfish and 
between shellfish, mud and vegetation percent 
reflectances. Although, the spectral curves of 
mud and shellfish are more alike than between 
vegetation and shellfish-mud spectral 
signatures. 

Comparing the individual spectral signatures 
of the shellfish shows a percent reflectance 
difference within the shellfish group. BOB4-1 5 
and 16 represent the “G” strata that are semi- 
wet with a thin film of mud. BOB4-1 and 2 
represent a stratum that is more like the “D” 
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strata with few vertical oysters and, some 
substrate showing and washed (dead) oysters 
on top of the shellfish reef that is dry. BOB4-1 
and 2 show a greater amount of reflectance 
that may be due in part to being drier and the 
shells being in a reclined position. These initial 
spectral signatures show collectively is there 
are appreciable differences not only between 
mud and shellfish but there is a degree of 
variability within shellfish strata. More research 
will be done during the course of this research 
to ascertain the significance of wet versus dry 
between shellfish strata types. 

Recent imagery of Murrells Inlet, SC, was 
acquired with the Airborne Imaging 
Spectrometer for Applications (AISA) sensor in 
September 2003. The AISA sensor has 35 
bands in the visible and near-infrared portion 
of the spectrum (441 - 873 nm). AISA 
bandwidths are either 6 nm (441 - 500 nm) 
or 3nm ( >500 nm) wide (Table 1). 


Figure 3. BOB4 monitoring site. 


Imagery was acquired by the CALMIT 
Hyperspectral Aerial Monitoring Program in 
Lincoln, Nebraska. Spatial resolution of the 
imagery was 3 m x 3 m. Four flight lines were 
needed to cover the entire estuary and part of 
the uplands. Imagery was flown within 2 hrs of 
low tide on September 20, 2003. Conditions 
were cloud free. Figure 4 depicts 6 of the 8 
flight line data files covering Murrells Inlet. The 
two files not shown cover upland areas to the 
west of the estuary. 

Imagery has been processed to scaled 
percent reflectance and georegistered. We are 
awaiting radiance files and updated, more 
accurately georegistered imagery. IKONOS 
imagery was available for 1999 but was not 
available in 2003 for the Murrells Inlet location, 
under the atmospheric conditions specified. 

Biophysical data was collected in situ , 24 hrs 
after the imagery collection. On September 
21 st , within a window of 2 hrs on either side of 
low tide, 6 different locations within the estuary 


were sampled. The locations represented 
areas adjacent to the varying conditions 
surrounding the estuary (developed, 
undeveloped, and to the center of the estuary) 
(Figure 5). At each of the 6 different locations, 
5 sampling sites were chosen for biophysical 
sampling. The 5 sites formed an L-shape in 
order to capture variation within the Spartina 
altemiflora (Figure 5). A transect along the 
high berms of the low marsh consisted of 2 
sample points (A and B). A Three additional 
samples (C, D, and E) were taken along a 
transect perpendicular to the first transect. The 
purpose of the second transect was to capture 
the variation in Spartina height from the low 
marsh along the stream berms to the higher 
elevations of the high marsh (Jensen et al., 
1998). Sampling locations were approximately 
10 meters apart. Each location was marked by 
a GPS point, taken with a Trimble GPS. 

The in-situ biophysical measurements 
obtained included total above ground biomass 
and leaf area index (LAI). Total above ground 
biomass was harvested for 0.25 x 0.25 m 
quadrats and placed into labeled plastic bags. 
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Table 1. Comparison of 
ATLAS Multi 




AISA Hyperspectral 


Spatial 

Resolution 

3 x 3 m 


Maximum 
number of 
channels 

35 


Band Number 

Band Center 

Bandwidth 

1 

441 

6 

2 

463 

6 

3 

474 

6 

4 

483 

6 

5 

500 

6 

6 

542 

3 

7 

550 

3 

8 

561 

3 

9 

585 

3 

10 

597 

3 

11 

609 

3 

12 

622 

3 

13 

632 

3 

14 

639 

3 

15 


3 

16 

670 

3 

17 

685 

3 

18 

692 

3 

19 

700 

3 

20 

705 

3 

21 

709 

3 

22 

714 

3 

23 

719 

3 

24 

728 

3 

25 

749 

3 

26 

757 

3 

27 

764 

3 

28 

774 

3 

29 

791 

3 

30 

808 

3 

31 

814 

3 

32 

819 

3 

33 

826 

3 

34 

849 

3 

35 

873 

3 
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the AISA Hyperspectral and 
spectral Sensors 



ATLAS Multispectral 


Spatial 

Resolution 

3.08 x 3.08 

meters 




9 



Band Number 

Band Center 

Bandwidth 














1 

485 

70 
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560 

80 
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615 
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Figure 4. The Airborne Imaging Spectroradiometer for Applications (AISA) sensor was used to 
acquire imagery of Murrells Inlet, SC on September 20, 2003. The 6 datasets covering the estuary 
are presented above. The 2 remaining datasets cover additional upland areas to the west of the 
estuary. Each image is presented in units of scaled percent reflectance in band combination 34, 18, 
2. Spatial resolution is 3 x 3m. North is oriented to the top of the page. 
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I 

I 


Figure 5. Ground referencing sampling pattern for the AISA hyperspectral data collection, 
September 20, 2003. The above sampling point location (site #7) representative of two transects for 
each of the six sampling locations: 1 ) a transect to capture the greater amounts of biomass along the 
berm, in the low marsh and 2) another transect to capture the transition in Spartina biomass from the 
low to high marsh, perpendicular to the berm. Inset graphic displays GPS points on an ADRG base 
map. AISA hyperspectral imagery displayed in band combination 32, 18, 2. 


I 
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LAI was measured with a Decagon 
Ceptometer and included readings above and 
below the Spartina altemiflora canopy (Jensen 
et al., 2002). Spartina spectra readings were 
obtained using a handheld GER 
spectroradiometer with a spectral range of 290 
to 1098 nm. Initial assessment of the LAI 
readings indicates capture of variation in 
Spartina LAI measurements from the high 
berms of the low marsh to the flat, slightly 
higher elevations of the low marsh (Figure 6). 
Harvested quadrats of Spartina were dried in a 
storage shed for approximately three months. 
The contents of each marked bag were 
transferred to another bag and weighed before 
being placed in a drying oven. Final dry weight 
was determined after 1 0 days. 

Statistical relationships between in situ and 
AISA-derived spectra, LAI, and biomass are 
being established. Using these mathematical 
relationships, vegetation indices based on 
previously conducted BRDF signatures of 
Spartina altemiflora will be used to generate a 
biomass map of the peak 2003 Spartina in 
Murrells Inlet (Jensen et al., 2002; Meisburger, 
2000; Shew and Linthurst, 1981; Thenkabail et 
al., 2000). 

The AISA spectra also will be convolved to 
mimic the ATLAS multispectral bandwidths. 
Using the previously derived vegetation 
indices applied to the ATLAS 1997 data, the 
convolved AISA 2003 data will be used to 
generate a biomass map. These two biomass 
maps will be compared. Of particular interest 
will be whether or not change in the Spartina 
may be correlated to changes or stress in the 
upland land cover (Porter et al., 1996; 
Vemberg et al., 1992). 

Further research questions will be addressed 
using this dataset. The AISA 2003 imagery 
and the convolved AISA imagery will be 
compared to quantify any improvement in the 
capture of Spartina altemiflora biophysical 
parameters from hyperspectral versus 
multispectral imagery. 
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Site # 4 
Site # 5 
Site # 6 
Site # 7 
Site # 8 
Site# 10 


Figure 6. When compared by quadrat, the trends in LAI are similar. The observed increase in 
LAI reading for quadrat sites “E” may indicate that the sampling scheme placed quadrats far enough 
into the marsh that they perhaps represent the topographic trend back to low” marsh near another 
creek. There was no sample site #s 1-3, 9. 
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Prentiss Hall, Saddle River, NJ, 576 pp. 


Published Papers 

Allen, J. and K. Lu (2003) Modeling and Prediction of Future Urban Growth in the Charleston Region 
of South Carolina: A GIS-based integrated approach. Conservation Ecology 8(2): 2. 

Allen, J., Kang S. Lu, and Thomas D. Potts. (2002) A GIS-based analysis and prediction of land-use 
change in a coastal tourism destination area. In Miller, Marc I., Auyong, Jan, and Hadley, Nina P. 
(eds.). Proceedings of the 1999 International Symposium on Coastal and Marine Tourism: 
Balancing Tourism and Conservation. (26-29 April 1999, Vancouver, BC, Canada). Washington 
Sea Grant Program and School of Marine Affairs, University of Washington, Seattle, WA, Oregon 
Sea Grant College Program, Oregon State University, Corvallis, OR, and Oceans Blue 
Foundation, Vancouver, BC. 

Aisup D.L., C.R. Coombs (1998) A Study of Plant Distribution Changes in Abandoned Rice 
Fields within Freshwater Estuaries, Georgetown-Horry Counties, South Carolina. 

Proceeding First International Conference on Geospatial Information in Agriculture and 
Forestry, pp. 11-653-11-659. 

Burdett L„ C. Runyon, W. McFee (2003) Mapping of Marine Mammal Entanglement Wounds: 
Sources of Mortality in Commercial Fisheries Predicted with GIS. ArcNews, Fall 2003. 
http://www.esri.com/news/arcnews/fall03articles/sources-of-mortalitv.html 

Bom, K., D.E. Porter, M. Neet and S. Walker. In press. An Internet mapping system for coastal 
resource management and data dissemination in support of a Special Area Management Plan. 
California and the World Ocean '02 Proceedings. 

Gardner, R. and D.E. Porter. (2001) Stratigraphy and geologic history of a southeastern salt marsh 
basin, North Inlet, South Carolina, USA. Wetlands Ecology and Management. 

Hulin A., C.R. Coombs, D.E. Porter (1999) A Policy Analysis of the Safe Harbor Act to Protect 
the Endangered Red-Cockaided Woodpecker Using GIS and Remote Sensing. Jour. 

Forested Wetlands. 

Jensen J.R., C.R. Coombs, D.E. Porter, B. Jones, S. Schill (1999) Smooth Cordgrass ( Spartina - 
altemiflora) Biomass and Leaf Area Index Parameters from High Resolution Imagery. 

Geocarto, v. 13, no. 4, pp. 25-34. 

Kleppel, G.S., D.E. Porter, M.R. DeVoe and R.H. Becker. Preserving biodiversity in rapidly urbanizing 
estuarine ecosystems. Science. 

Nichols, C.C., T.C. Siewicki, and J.W. Daugomah, D.E. Porter, B. Jones and J.R. Jensen. 2000. 
Coastal land-cover classification using NASA’s Airborne Terrestrial Applications Sensor (ATLAS) 
data. Proceedings of Sixth International Conference on Remote Sensing for Marine and Coastal 
Environments, pp. 279-286. 
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Savitsky B.G., J.S. Allen, K.F. Backman (1999) The role of Geographic Information Systems in 
Tourism Planning and Rural Economic Development. Tourism Analysis. Sagamore 
Publishing Co., Champaign-Urbana, IL 

Schill S.R., J.R. Jensen, D.E. Porter, B.C. Jones, and C. Coombs (2000) A discriminant 
analysis on the biophysical characteristics of Spartina altemiflora and muitispectral 
reflectance in remotely sensed data. Proc. 6* International Conference on Remote 
Sensing for Marine and Coastal Environments. May 1-3, Charleston, SC. pp. 89-102. 

Schill, S.R., J.R. Jensen, G.T Raber, D.E. Porter, 2004, Temporal Modeling of Bidirectional Reflection 
Distribution Function (BRDF) in Coastal Vegetation, GIScience and Remote Sensing. 41 (1 ):61- 
80. 

White, D.L., D. Bushek, D.E. Porter and D. Edwards. (1999) Geographic Information Systems 
(GIS) and Kriging: Analysis of the Spatial and Temporal Distributions of Perkinsus marinus in 
a Developed and an Undeveloped Estuary. Journal of Shellfish Research. V. 1 7, no. 5, pp. 
1473-1476. 

White, D., D.E. Porter and A.J. Lewitus. 2004. Spatial and temporal analyses of water quality and 
phytoplankton biomass in an urbanized versus pristine salt marsh estuary. Journal of 
Experimental Biology and Ecology. 298:255-273. 

Vincent, J.S., D.E. Porter, L.D. Coen, D. Bushek and S. Schill. 2003. Using hyperspectral remote 
sensing to map and assess intertidal shellfish resources in the southeastern US. Journal of 
Shellfish Research. Vol. 2(1 ):359. Reviewed abstract 


Presentations 

Conference Presentations (abstracts submitted for each) 

Allen, J. (2000) Land-Use Change and Urban Growth Modeling for Coastal South Carolina.” South 
Carolina State Mapping Advisory Committee Biennial Conference. Columbia, South Carolina. 
January 26-27, 2000. 

Allen, J. and Lu, K. S. (1999) A GIS-Based Analysis and Prediction of Land-Use Change in a Coastal 
Tourism Destination Area. 1999 World Congress on Coastal and Marine Tourism. Vancouver, 
British Columbia, Canada. April 26-29, 1999. 

Allen, J.S. and S.L. Kang (1998) "Population and tourism growth and its relation to ecosystem 
health" Conference on Southeast Coastal Ocean Research. Savannah, GA., April 1998. 

Alsup D.L., C.R. Coombs, M. Weiner (1998) The monitoring of wetland ecosystem health in 
impounded estuary communities using remote sensing. College of Charieston, April 1997. 

Bom, K„ D.E. Porter, M. Neet and S. Walker (2003) An Internet mapping system for coastal resource 
management and information dissemination in support of a Special Area Management Plan. 
Coastal GeoTooIs '03. Charleston, SC. January 2003. 

Burdett.L., W. McFee, C. Runyon (2004) Working Backwards: Using the rope impressions of 

entangled marine mammals to create GIS maps that identify potential sources of entanglement 
Southeast and Mid-Atlantic Marine Mammal Symposium held March 26-28 at Harbor Branch 
Oceanographic Institution, N.C. 
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Bushek D., D. White, D.E. Porter, B. Jones, J. Keesee and D. Edwards (1998) Relationships 
between land-use patterns and dermo disease in oysters from two South Carolina 
estuaries. 3 rd International Symposium on Aquatic Animal Health. Baltimore, MD. August, 

1998. 

Bushek, D., D. White, D. Porter and D. Edwards (1998) "Land-use patterns, hydrodynamics and 
the spatial pattern of Dermo disease in two South Carolina estuaries’. Aquaculture '98, Las 
Vegas, NV, February 1998. 

Chose J.R. (1998) Tidal Creek Bank Erosion in Murrells Inlet, SC. Southeastern ARC/INFO 
Users Group Meeting, Charleston, SC, October, 1998. 

Chose J.R. and C.R. Coombs (1998) Parameters influencing creek bank erosion in two South 
Carolina tidal salt marshes. S.C. Acad. Sciences, Clemson, SC. March 1998. 

Chose J.R. and C.R. Coombs (1998) Creek bank erosion of the tidal salt marshes of Murrells 
Inlet and North Inlet, South Carolina. Conference on Southeast Coastal Ocean Research. 
Savannah, GA., April 1998. 

Coombs C.R., J.R. Chose, D.L. Alsup, and J. Daugomah (1997) The monitoring of ecosystem health 
in South Carolina estuaries using remote sensing: A NASA EPSCoR Project. College of 
Charleston, April 1997. 

Daugomah J. (1998) Natural and Anthropogenic Factors Affecting Spatial Palaeomonetes pugio 
Population Abundance in Murrells and North Inlets, SC. Southeastern ARC/INFO Users Group 
Meeting, Charleston, SC, October, 1998. 

Goddard M. (2000) Integrating Landsat Thematic Mapper Data and BASINS Integrated Modeling for 
Water Quality Mapping and Prediction. ASPRS Regional Meeting Student Presentation, 

Clemson University. 

Hay G. (2001) Modeling the Habitat Characteristics of Tidal Creeks Used by Recreationally Important 
Juvenile Finfish with a Geographic Information System. Coastal GeoTools 2001. Charleston, SC. 
January 2001. 

Hulin A. (1998) Endangered Red-Cockaided Woodpecker’s Foraging and Nesting Habitat. 
Southeastern ARC/INFO Users Group Meeting, Charleston, SC, October, 1998. 

Kleppel G., D.E. Porter, J. Allen and M.R. DeVoe (1999) Implications of changing demographic and 
land use patterns to coastal ecosystem integrity in the southeastern United States. 9 th Annual 
Meeting of SET AC-Europe. Leipzig, Germany. May, 1999. 

Meisburger, J.L., J.R. Jensen, S. Schill and D.E. Porter. (2000) Quantification of biomass and leaf are 
index in a Charleston, SC estuary using low-altitude AVIRIS data. American Society for 
Photogrammetry and Remote Sensing Annual Convention. Washington, DC. May 2000. 

Nichols C.C., T.C. Siewicki, and J.W. Daugomah, D.E. Porter, B. Jones and J.R. Jensen. (2000) 
Coastal land-cover classification using NASA’s Airborne Terrestrial Applications Sensor (ATLAS) 
data. Sixth International Conference on Remote Sensing for Marine and Coastal Environments. 
Charleston, SC. May 2000. 

Olson, G., J. Jensen, S. Schill and D.E. Porter. (2001) Remote sensing of wetland biophysical 
variables in the ACE Basin National Estuarine Research Reserve. Coastal GeoTools 2001. 
Charleston, SC. January 2001. 
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Porter, D.E. (2000) Using innovative remote sensing and image processing techniques for assessing 
coastal ecosystem health. NASA/EPSCoR Earth Science Enterprise Workshop. Clemson, SC. 
August 2000. Invited presentation. 

Porter D.E., G.T. Chandler, B.E. Jones and M.W. Williamson (1998) The use of spatial modeling to 
map estuarine areas susceptible to photo-induced toxicity. Southeastern ARC/INFO Users 
Group Meeting, Charleston, SC, October, 1998. 

Porter, D.E. (1999) An overview of the role of geographic information processing techniques in 

support of coastal environmental decision making. Invited presentation. Southeast Partnership for 
Environmental Technology Education Fourth Annual Resource Conference. Charleston, SC. April 
1999. Invited presentation. 

Porter D.E., D. Bushek, D.L. White and J. Keesee (1998) Host-Parasite Relationships and Disease 
Status as a Measure of Ecosystem Health. Conference on Southeast Coastal Ocean Research. 
Savannah, GA., April 1998. 

Porter, D.E., J. Allen, T. Siewicki, M. Gielazyn, D. Edwards and W.K. Michener (1998) Defining the 
Role of Database Management in Support of Long-term, Multidisciplinary Environmental 
Research Efforts. Conference on Southeast Coastal Ocean Research. Savannah, GA., April 
1998. 

Runyon, C.J., L. Burdett, W. McFee (2003) Using GIS as a Forensic Tool to Predict Sources of 
Marine Mammal Entanglement in Commercial Fisheries. The 15th Biennial Conference on the 
Biology of Marine Mammals, Joseph S. Koury Convention Center, Greensboro, NC. December 
14-19, 2003. 

Schill, S., J. Jensen and D.E. Porter. 2001. Measurement and modeling of the spectral and directional 
reflection properties of smooth cordgrass ( Spartina altemiflora). Coastal GeoTools 2001. 
Charleston, SC. January 2001 . 

Schill, S., J. Jensen and D.E. Porter. 2001, Temporal modeling of BRDF patterns and phenological 
change in coastal vegetation. 2001 American Society of Photogrammetry and Remote Sensing 
Annual Convention and Exposition. St. Louis, MO. April 2001. 

Schill, S.R., D.E. Porter and J.R. Jensen. 2000. Acquiring and modeling BRDF in coastal 
environments. Mid-South Region of the American Society of Photogrammetry and Remote 
Sensing Annual Meeting. Clemson, SC. October 2000. 

Schill S.R., J.R. Jensen, D.E. Porter, B.C. Jones, and C. Coombs. 2000. A discriminate analysis on 
the biophysical characteristics of Spartina altemiflora and multispectral reflectance in remotely 
sensed data. Sixth International Conference on Remote Sensing for Marine and Coastal 
Environments. Charleston, SC. May 2000. 

Schill, S.R., D.E. Porter, and J.R. Jensen. 2000. Evaluating bi-directional reflectance distribution 
function (BRDF) of smooth cordgrass ( Spartina altemiflora). National Estuarine Research 
Reserve Association Annual Meeting. Williamsburg, VA. October 2000. 

Vaughan, D.F. and J. R. Jensen. Extraction of Spartina altemiflora biophysical parameters using 
narrow-band vegetation indices and high spatial resolution hyperspectral imagery. Association of 
American Geographers 2004 Annual Conference, March 14-19, 2004, Philadelphia, 

Pennsylvania. 

Vaughan, D.F., J.R. Jensen, D.E. Porter, S. Walker, B. Hadley, and J. Tullis. Hyperspectral Change 
Detection of Spartina altemiflora in Murells Inlet, South Carolina. American Society of 
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Photogrammetry and Remote Sensing 2004 Annual Conference, May 23-28, 2004, Denver, 
Colorado. 

Vemberg, F.J., W.B. Vemberg, D.E. Porter, G.T. Chandler, H.N. McKellar, D. Tufford, T. Siewicki, M. 
Fulton, G. Scott, D. Bushek and M. Wahl. Impact of coastal development on land-coastal waters. 
1999. Medcoast 99. Antalya, Turkey. November 1999. 

Vincent, J.S., D.E. Porter, L. Coen, D. Bushek and S. Schill. 2003. Assessing and mapping shellfish 
resources using hyperspectral remote sensing. Coastal GeoTools ’03. Charleston, SC. January 
2003. 

White, D.L., D.E. Porter and A.J. Lewitus. (2001) Spatial heterogeneity of nutrients and primary 
production in an urbanized estuarine system. Coastal GeoTools 2001. Charleston, SC. January 
2001. 

White, D., D.E. Porter, and A.J. Lewitus. (2000) Spatial and temporal analyses of water quality 
patterns in two well-mixed, high-salinity estuaries. National Estuarine Research Resen/e 
Association Annual Meeting. Williamsburg, VA. October 2000. 

White, D.L., D.E. Porter, A.J. Lewitus, M. Neet, D. Tufford and J.R. Jensen. (2000) Integrating in situ 
sampling and remote sensing techniques to assess and model spatial heterogeneity of nutrients 
and primary production in an urbanized estuarine system. Mid-South Region of the American 
Society of Photogrammetry and Remote Sensing Annual Meeting. Clemson, SC. October 2000. 

White, D., D.E. Porter, J. Jensen, S. Schill and B. Jones. 1999. An evaluation of non-destructive 
techniques for estimating aboveground biomass of the cordgrass Spartina altemiflora. Iff* 
International Biennial ERF Conference. New Orleans, LA. September 1999. 

White D., D.E. Porter and D. Bushek (1997) A Comparison of Spatial and Temporal Patterns of 
Perkinsus marinus in Oyster Populations in a Developed Estuary and an Undeveloped Estuary in 
South Carolina. Spatial Data and Remote Sensing in Invertebrate Fisheries Habitat, Research 
and Management Workshop. Fort Walton Beach, FL. 19-20 April 1997. 


Invited Presentations 

Allen J.S. (2001 ) Growth Modeling and Growth Management in South Carolina. Interview with News 
Channel 7, Spartanburg, SC. January 2001. 


Allen J.S. (2000) Using a Growth Prediction Model for Monitoring Coastal Ecosystem Health. 
NASA/SC EPSCoR project workshop. Clemson, South Carolina. August 21 , 2000. 

Allen J.S. (2000) Growth Management - What it Means to the Palmetto State. South Carolina 

Association of Conservation DistrictsAJSDA-NRCS/SCDNR Joint Conference. Charleston, South 
Carolina. January 21, 2000. 

Allen J.S. (1999) Growth Prediction Modeling: Lessons for Upstate South Carolina. Reedy River 
Growth Management Conference. Greenville, South Carolina. December 1, 1999. 

Allen J.S. (1999) Forces Behind Land Use Planning: A GIS-Based Growth Prediction Model.” 31st 
Annual Meeting of the South Carolina Forestry Association. Myrtle Beach, South Carolina. 
November 3-5, 1 999. 
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Allen J.S. (1999) Forces behind land use planning: A GIS-based growth Prediction model. Annual 
S.C. Forestry Association, Myrtle Beach, S.C., November, 1999. 

Allen J.S. (1 999) GIS and Predicting Growth. Annual Conference of the South Carolina Association of 
Special Purpose Districts, Myrtle Beach, S.C., October, 1999. 

Allen J. and S.L. Kang (1999) GIS-based Analysis and Prediction of Land-Use Change in Coastal 
Tourism Destination Areas. Proceedings of World Congress on Coastal and Marine Tourism. 
April 29, 1999. Vancouver, British Columbia, Canada. 

Coombs [Runyon] C.R. (1999) Monitoring Wetlands and Ecosystem health in coastal South Carolina 
using GIS and remote sensing, Denver, CO. Annual GSA meeting, October 1999. 

Coombs [Runyon] C.R. (1999) Modeling Coastal Ecosystem Health Using Advanced Techniques of 
Remote Sensing and GIS. Southeastern ARC/INFO Users Group Meeting, Charleston, SC, 
October, 1998. 

Coombs C.R., J.R. Jensen, D.E. Porter (1998) Use of remote sensing and GIS as monitoring tools for 
wetlands and ecosystem health. Texas Tech University, Lubbock, TX, February 20, 1998. 

Coombs C.R. (1997) A new way to explore wetlands in the classroom. KidSat (EarthCam) Training 
Workshop, College of Charleston, July, 1997. 

Jensen, J.R., D.E. Porter and C.R. Coombs (1998) Biophysical Remote Sensing of Coastal 

Wetlands. AMS Annual Meeting and Science Innovation Exposition. Philadelphia, PA. February 
1998. 

Porter D.E. (1999) An overview of the role of geographic information processing techniques in 
support of coastal environment decision making. Southeast Partnership for Environmental 
Technology Education Fourth Annual Resource Conference. Charleston, SC, April 1999. 

Porter, D.E. (1998) Murrells Inlet Under the Microscope: Preliminary Findings and Recommendations 
for Coastal Zone Management. Murrells Inlet 2007 Task Force and community of Murrells Inlet. 
January 1998. 


Award-winning presentations 

Allen J., Lu K.S. (2000) Most innovative use of technology in S.C. Awarded by SC Chapter of the 
American Planning Association for Charleston Area Urban Growth Model. 

Dr. J. R. Jensen, 2004, John E. Estes Award, ASPRS, Sponsored by SAIC. American Society of 
Photogrammetry and Remote Sensing, Denver, Colorado 

Porter D.E., G.T. Chandler, B.E. Jones and M.W. Williamson (1998) The use of spatial modeling to 
map estuarine areas susceptible to photo-induced toxicity. Southeastern ARC/INFO Users 
Group Meeting, Charleston, SC, October, 1998. 

Schill S. (2000) BRDF Modeling in estuarine systems. Remote Sensing Division, ASPRS 2000, 
Washington D.C. 1 st Place. 

White D. Geographic Information Systems (GIS) and Kriging: Analysis of the Spatial and Temporal 
Distributions of the Oyster Pathogen Perkinsus marinus in a Developed and an Undeveloped 
Estuary. Southeastern ARC/INFO Users Group Meeting, Charleston, SC, October, 1998. 
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Section Ilia: Spin-Off Projects (Funded) 


Funding Period 
Amount ($) 

Granting Agency 


Project Title 

2004-2009 

$500,000 

SC DNR 

J. Allen et al. 

Reedy River Watershed Education Program 

2003-2005 

VKR Foundation 

J. Allen et al. 

The Saluda-Reedy Watershed Consortium: 

$270,000 



Taking Action for Water Quality Improvement 
and Watershed Management 

2003-2007 

$2.4 million 

NASA/REASoN 

J. R. Jensen 
M.E. Hodgson 
S. Cutter 
D. J. Cowen 

Development of Remote Sensing-assisted 
Natural and Technological Hazards Decision 
Support Systems 

2004-2005 

$158,140 

NOAA CICEET 

D. E. Porter 
J.R. Jensen 
V. Klemas 

NERRS Remote Sensing Applications 
Assessment Project (RESSAP)". 

2002-2004 

$245,491 

NOAA CICEET 

S. Schill 
D. E. Porter 
L. Coen 

Development of an Automated Mapping 
Technique for Monitoring and Managing 
Shellfish Distributions 

2000-2001 

$113,130 

NOAA 

Porter et al. 

Urbanization and Southeastern Estuarine 
Systems 

2000-2004 

NOAA Coastal 
Ocean Program 

Porter 

Development of a GIS-based Database 
Management Program to Characterize 
Sources and Effects of Natural Parameters 

$617,359 



and Anthropogenic Impacts to Coastal 
Ecosystems 

2000-2001 

$19,500 

NSF 

Porter et al. 

Multidimensional Consequences of Urban 
Encroachment on Natural Ecosystems 

Biid 

NASA/EPSCoR 

Porter 

Graduate Student Training Grant 

1999-2002 

$48,572 

EPA 

Porter et al. 

CISNet: Molecular to Landscape-scale 
Monitoring of Estuarine Eutrophication 

1998 - 1999 

$125,000 

NOAA Coastal 
Ocean Program 

Porter et al. 

Geographic Information Processing 
and Spatial Modeling on Urbanization 
in Southeastern Estuarine Systems 

1997 - 1998 

SC Sea Grant 

Porter, Allen 

A Spatial Data Assessment for LU- 

$33,744 



CES 


Section lllb: Spin-Off Projects (Proposed / Pending) 


Funding Period/ 
Amount 


PKs) 

Project Title 

2004 

$2,975,576 

NSF 

J. Allen et al. 

Biocomplexity as an Integrative Concept for 
Interdisciplinary Graduate Education and 
Research 

2004 

$410,922 

NSF 

J. Allen et al. 

impacts of Urbanization on social, economic 
and environmental values - Does typology 
matter? 

2004 

$2,000,000 

NSF 

J. Allen et al. 

Multiple Information Streams for Better 
informed Policy Makers 
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Section IV: Participants 


University of South Carolina (J. Jensen, D. Porter) 


Student 

EsnnsHi 

1 »J B 1 IlfWHBB 

Thesis/Dissertation Title 

George T. Raber 

Ph.D., 2003 

Geography 

The Impact of Varied Nominal Posting Density 
Lidar Data on DEM Accuracy, hydraulic modeling 
and Flood Zone Delineation 

Jason A. Tullis 

Ph.D., 2003 

Geography 

Data Mining to Identify Optimal Spatial 
Aggregation Scales and Input Features: Digital 
Image Classification with Topographic Lidar and 
Lidar Intensity Returns 

Steve R. Schill 

Ph.D., 2001 

Geography 

Modeling hyperspectral Bidirectional Reflection 
Distribution Function (BDRF) data of smooth 
cordgrass (Spartina altemiflora) using a 
sandmeter field goniometer 

Jeff Vincent 

Ph.D. 

Candidate 

Geography 

Mapping Shellfish Distribution Using 
Hyperspectral Remote Sensing 

David F. Vaughan 

Ph.D. 

Candidate 

Geography 

Examination of Hyperspectral Data for 
Classification and Biophysical Parameter 
Measurement of Smooth Cordgrass (Spartina 
altemiflora) 

Brian Gray 

Ph.D. 

Candidate 

Biostatistics & 
Epidemiology 

Modeling spatially and temporally correlated 
oyster infection data 

Mark Jackson 

Ph.D. 

Candidate 

Geography 

On leave of absence 

Jennifer Meisberger 

Ph.D. 

Geography 

Quantification of Biomass and Leaf Area Index 
(LAI) Using AVIRIS imagery 

David L. White 

Ph.D. 

Candidate 

Marine Science 

Land Cover and Water Quality: A Spatial 
Investigation of Phytoplankton community 
structure in two well mixed marsh ecosystems 

Xinghe Yang 

Ph.D. 

Geography 

Georeferencing CAMS Data: Polynomial 
Rectification and Beyond 

Yingming Zhou 

Ph.D. 

■eUJ.y.I.M 

Modeling Uncertainty of Areai Interpolation 

Katherine Bom 

M.S. 

Marine Science 

Using GIS, Remote Sensing and Internet 
Mapping Systems for Mapping, Assessing and 
Monitoring Localized Coastal Estuaries 

Gunner Olson 

M.S. 

Geography 

Remote Sensing of Estuarine Salt Marsh: 
Assessment of Biophysical Variables Utilizing 
High Resolution Data 

Judith Burgland 

M.S. 

Geography 

Evaluating the Use of Radarsat SAR and JERS-1- 
SAR for Estimating LAI and Delineating Upland 
Land Cover and Tidal Freshwater Wetland 
Vegetation Species in the Tivoli North Bay Area, 
New York 

Laura Schmidt 

M.S. 

Geography 

Evaluation of the Utility and Accuracy of LIDAR 
and IFSAR Derived DEMs for Flood Plain 
Mapping 

Jason Tullis 

M.S. 

Geography 

Artificial Intelligence for Automated Housing 
Enumeration from Pan-Sharpened IKONOS 
Imaqery 

George Raber 

M.S. 

Geography 

Vegetation Characterization and Improvement of 
Digital Surface model (DSM) Creation Using 
Multiple Return LIDAR 


I 

I 
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Student 

Degree 

E0S3333BBBI 

Thesis/Dissertation Title 

Matthew Neet 

M.S.P.H. 

Environmental 
Health Sciences 

Environmental Health Sciences Prototype 
Development of a GIS-based Web Presentation 
for Spatial Analysis and Promoting Increased 
Awareness of Fish Consumption Advisories for 
Lake Hartwell, SC 

Margaret 
(Williamson) Vo 

M.S.P.H. 

Environmental 
Health Sciences 

Integrating meiobenthic copepod bioassays and 
spatial modeling techniques to examine 
photo-enhanced toxicological hazard of MurTells 
Inlet sediments contaminated with polycyclic 
aromatic hydrocarbons 

Heath Kelsey 

M.S.P.H. 

Environmental 
Health Sciences 

Examination of Land Use and Fecal Coliform 
Bacterial Pollution From Human and Animal 
Sources in Murrells Inlet, South Carolina 

Cameron Kerr 

M.S.P.H. 

Environmental 
Health Sciences 

Using LIDAR in the development of watersheds 
and sub-watersheds for assessment and 
modeling of source loadings in natural resource 
management practices 

Giannina Dimaio 

M.S.P.H. 

Environmental 
Health Sciences 

A Review and Analysis of SC DHEC's OCRM 
Regulations for Permitting in Critical and Non- 
critical areas 

Katherine Bom 

M.S.P.H. 

Environmental 
Health Sciences 

Development of a GIS-based resource 
management tool to assess anthropogenic 
activities and physiographic influences on the 
health and structure of the Murrells Inlet, South 
Carolina estuary 





Cameron Kerr 

as. 

Marine Science 

N/A 

Giannina DiMaeo 

B.S. 

Marine Science 

N/A 


B.S. 

Marine Science 

N/A 
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Clemson University (J. Allen 


Student 


Thesis/Dissertation/Internship Title 


Kang Shou Lu 

Ph.D. 

Candidate 

Jeffrey Allen 

Ph.D. 

Candidate 

Megan Goddard 

M.S. 

Heather Crean 

M.S. 


Parks, Recreation ^acro-level Land-Use Change Analysis and 
& Tourism Mamt Pre diction for Coastal Tounsm Destination 
9 Areas 


Policy Studies 


Urban Growth Modeling and its Policy 
Implications 


Land Use Changes and its affects on water 
quality in a coastal watershed 


D1 ; e . .. Coastal South Carolina’s Population and 

Planning Studies Deve|opment Trends 


College of Charleston (C. Runvon f Coombs 


Student 


Thesis/Dissertation/Internship Title 


Marian Mitchell 

MES 

(2003) 

Policy / Geology 

Using Geographic Information Systems Software 
to Analyze the Capacity of Wood Stork Nesting 
Colonies in Relation to Wetland Availability 

Joy Parikh 

MES 

(2003) 

Geology 

Integrating AutoCAD Facility Data into an 
Environmental GIS Database for Impact 
Assessment Using Marshall Space Flight Center 
as a Model 

Robin Adams 

MES 

(2002) 

Policy / Geology 

Modeling and Restoration of Water Quality in 
Two Tidal Ponds on Kiawah Island, South 
Carolina 

Connie Moy 

(MES) 

Marine Biology 

Development and Validation of an Estuarine Biotic 
Integrity Index for South Carolina Tidal Creeks 

Gretchen Hay 

MES 

(2001) 

Geology 

Modeling the Habitat Characteristics of Tidal 
Creeks Used by Recreationally Important 
Juvenile Finfish with a Geographic Information 
System 

Scott Rutzmoser 

MES 

(2001) 

Geology 

Assessing the percent cover of Phragmites 
austalis in the Ace Basin NERR, SC 

Andrew Hulin 

MES 

(2000) 

Policy/Geology 

Evaluating the Use of GIS and Remote Sensing 
for Use in Establishing Land-Use for Endangered 
Species Habitat 

James Daugomah 

MES 

(2000) 

Geology 

Natural and Anthropogenic Factors Affecting 
Spatial Palaeomonetes pugio Population 
Abundance in Murrells and North Inlet, SC 

Dianna Alsup 

MES (1999) 

Geology 

Succession in Abandoned Rice Fields: A Study of 
Changes to Plant Distribution in a Marsh 
Environment 

Jaime Chose 

MES (1999) 


Parameters Influencing Bank Erosion in Two South 
Carolina Tidal Salt Marshes 


Scott Bruemmer B.S. 


Prentiss Lund B.S. 


Siobhan O'Reilly- B.S. 

Green Candidate 
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Section V: Collaborations 


NASA Centers 

NASA Goddard Space Flight Center 
NASA Stennis Space Flight Center 
NASA Johnson Space Flight Center 
NASA Langley Space & Research Center 
NASA Marshall Space Flight Center 

Business and Industry 

Calmit Hyperspectral Aerial Monitoring Program 
The Nature Conservancy 

SMARTech, Charleston, SC - satellite ground station, Small Business 
Geometries, Myrtle Beach and Columbia, SC - geotechnical/civil engineering, Small 
Business 

Other Federal 

National Oceanic and Atmospheric Administration 
Corps of Engineers 

State Agencies 

South Carolina Department of Natural Resources 

South Carolina Department of Health and Environmental Control 

South Carolina Department of Commerce 


Non-government Organizations 

Murrells Inlet 2007 Task Force 
Charleston City Council 


South Carolina EPSCoR - Wetlands Research Cluster 
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Section VI: Additional Uses of Data derived from this project 


Classes using data 

University of South Carolina (Jensen/Porter) 

Remote Sensing 751 
Remote Sensing 551 

Resource Management and Decision Making ENHS 775 

College of Charleston (Runyon) 

Geomorphology GEOL225 
Introduction to Remote Sensing GEOL314 
Advanced Remote Sensing GEOL442 
Geographic Information Systems GEOL449 

University of Charleston (Runyon) 

Applications in Remote Sensing EVSS642 
Geologic Applications in GIS EVSS 649 


South Carolina EPSCoR - Wetlands Research Cluster 
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Appendix 1: 


Land Use Modeling and Prediction Using an Artificial Neural Network 


Jeffery S. Allen and Kang S. Lu 

Strom Thurmond Institute, Clemson University, Perimeter Road, Clemson, SC 29634; 
email: ieff@strom.clemson.edu klu@strom.demson.edu 


1 Introduction 

One of the challenges of urban growth simulation is to identify and determine the rules to be used 
in a model that fit the unique geographic environment (Batty, 1998; Wu, 2002). Regardless of the 
approach and the form that a model may take, these rules are either conceptually preset by urban 
experts or empirically derived from case studies. Experts such as Von Thunen (1826) in the field 
of urban economics strived to use deterministic rules and a select few variables in modeling 
urban space. Their models succeeded to varying degrees in explaining some spatial phenomena 
but failed in practical applications because they all oversimplified the urban reality (Lee, 1973; 

Sui, 1997). The complexity of urban systems also makes it difficult, if not impossible, to derive 
universal rules for building a knowledge-based expert system that is applicable everywhere like 
those in the medical science even though some attempts have been made in this direction 
(Pijanowski et al., 1997, Allen and Lu, 2003). It is only possible or more appropriate to empirically 
derive the regionalized or localized rules that govern each of the spatially differentiated urban 
systems. 

All empirical models rely on sample data for training or calibrating and assume that rules or 
relationships, often in terms of coefficients or weights, revealed from these data remain to be true 
within a certain geographic domain for a given time period. Although there are many methods 
available for deriving the coefficients or weights for conceptual models, they mainly fall into two 
groups: conventional statistical techniques and innovative neural networks. Among statistical 
models, logistic regression has been popular for several reasons (see McMillan, 1989; Landis & 
John, 1997; Allen, Lu, and Potts, 2002; Wu, 2002; Allen & Lu 2003). According to Allen and Lu 
(2003), the logistic framework is capable of handling discrete land use variables and the mix of 
both discrete and continuous independent variables; it is nonlinear in form to better represent the 
nature of the complex urban reality; it is flexible to be tailored to individual land use systems; and 
it is available within most statistical packages. However, due to both the complexity of urban land 
use systems and limitations of the model, it does not always provide satisfactory predictions, as 
several studies have shown (Landis, 1994, Allen & Lu, 2003). Like other mathematical models, 
the model assumes that all variables be independent from one another, which appears to 
contradict the hierarchical urban reality in which all factors are interrelated and interdependent. 
This leads to a search for alternative models such as neural networks based on the generic 
algorithms. 

This past year, the land use project team explored the applicability of artificial neural networks in 
land use simulation, developed an operational neural net model for predicting urban growth, 
applied the model in the Myrtle Beach region of South Carolina, and tested its reliability and 
validity using both areal sample datasets and global sample datasets against the logistic 
regression. 


2 Artificial neural networks 

Neural network computing and generic algorithms have been proposed as possible alternative 
approaches to handle the uncomputability of complex systems (Sui 1997). Although the first 
neural networks model appeared in the early 1940s (McCulloch & Pitts, 1943), they were mainly 
applied in signal processing (Widrow & Steams, 1985), control (Nguyen & Widrow, 1989, Miller, 
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Sutton, & Werbos, 1990), pattern recognition (Le Cun et at., 1990), medicine (Anderson, 1986; 
Anderson, Golden, and Murphy, 1986; Hecht-Nilsen, 1990), speech production and recognition 
(Sijowski and Rosenberg, 1986; Lippmann, 1989), and business (Collins, Ghosh, and Scofield, 

1 988). Like the history of land use modeling, neural network studies were few in a period of quiet 
years in the 1 970s after the first golden age of the 1 950s and 1 960s, but revitalized vigorously in 
the late 1980s and 1990s mainly due to the discovery of robust propagation training methods and 
the advent of computer technology (Fausett, 1994). However, it was not until the mid-1990s that 
neural networks were used in the field of land use modeling and resource management. 

Openshaw (1993) started using a neural net for modeling spatial interaction. Wang (1994) used 
artificial neural networks in a geographical information system for agricultural land-suitability 
assessment. Gimblett et al. (1994) developed a forest management decision model based on 
neural networks and generic algorithms and tested the model in the Hoosier National Forest. 
Fischer and Gopa (1994) used an artificial neural network approach to the modeling of 
interregional telecommunication flows. Gong (1996) combined evidential reasoning and artificial 
neural network techniques for geological mapping. More recently, Yeh and Li (2002) used neural 
networks and cellular automata to simulate potential or alternative urban development patterns 
based on different planning objectives. 

According to Kohonen (1988), neural networks are “massively parallel interconnected networks of 
simple (usually adaptive) elements and their hierarchical organizations which are intended to 
interact with other objects of the real world in the same way as the biological nervous systems 
do.” The neural network model has certain advantages over mathematical methods (linear and 
nonlinear) in that it can handle interdependency among the factors. The model also differs from 
conventional rule-based expert systems in that rules governing the connections among neurons 
are generated by learning from examples rather than by the preset knowledge from experts. 
Previous studies indicate that neural networks provide superior levels of performance to those of 
conventional statistical models because they can well handle the uncertainties of spatial data. 
Neural networks also are open and flexible compared to other deterministic models in that they 
can be tailored to the available variables and datasets. 

There are many types of neural networks that have been developed for addressing different 
problems. One of them is the backpropagation neural net that has the capacity to recognize and 
classify patterns through training or learning processes (Yeh and Li 2002). Figure 1 shows a 
multiplayer neural network. The net contains an input layer with multiple units, a hidden layer with 
multiple units, and an output layer with multiple units. Units between two adjacent layers are 
interconnected. Each unit from the input layer sends its signal to each unit of the hidden layer. 
Each hidden unit receives the signal from each input unit and sums the signal with different . 
weights and applies an activation function and sends its own signal to each output unit. Each 
output unit receives a signal from each hidden layer and sums the signals with corresponding 
weights and computes the output. This process can repeat if there are more hidden layers. The 
weights can be determined using the robust backpropagation algorithm to be discussed next. 
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Backpropagation Neural Network 
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Figure 1. A multi-layer backpropagation neural network. 


3 Backpropagation training algorithm 

The training of a network by backpropagation involves three stages: the feedforward of the input 
training pattern, the backpropagation of the associated error, and the adjustment of the weights. 
There are many methods available. The algorithm described below is for a standard three layer 
backpropagation neural net. 

The training starts with initialization of weights and bias by assigning a small random number to 
each of them. This can be achieved using the technique developed by Nguen-Widrow to keep 
random values between -0.5 and 0.5, as is commonly the case. For each record or training 
pattern in the dataset, the training will complete the following three stages. This will continue in an 
iterative fashion until the stop condition is met. 


3.1 Feedforward training 

During the feedforward, each input unit (X/, /= 1, ..., n) receives signal x, and broadcasts this 

signal to each of the hidden units. Each hidden unit (Z , j = 1 p) sums its weighted input 

signals, 


/= i 
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z J n r v <* + 1 V- 

applies its activation function to compute its output signal, 

z f f ( z J n j), 

and sends this signal to all units in the output layer. Each output unit (Y k , k= “\, .... m) sums its 
weighted input signals, 

yjn k =w 0k + {w jk z j 
7=1 

applies its activation function to compute its output signal, 

y k =f(yjn k ) 

which is the response of the net for the given input pattern. 

3.2 Error Backpropagation 

During training, each output units compares its computed activation y k with its target value t k to 
determine the associated error for that pattern with that unit. Based on this error, the factor 6 k ( k 

1, .... m) is calculated. Each output unit (Y k , k= 1 m) receives a target pattern (t k ) 

corresponding to the input training pattern, computes its error information term 

S k =(t k -y k ) f'fyjnj, 

calculates its weight correction term (used to update w jk later) with a learning rate a, 

Aw jk = a 5 k Zj, 

calculates its bias correction term (used to update iv a later), 

& w ok = a z j< 

and sends to units in the layer below (the hidden units). 

Each hidden unit (Z-,/ = 1, .... p) sums its delta inputs (from units in the layer above) 

5J n r I 5 k w ,k ■ 

k= l 

multiplies by the derivative of its activation function to calculate its error inform term, 

fff 5 - in i f' (zJn,), 

calculates its weight correction term (used to update later) with a learning rate a , 

calculates its bias correction term (used to update w 0k later), 

%= a 5 r 
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3.3 Updating weights and bias 

Once the weight adjustment terms are calculated, they can be used to update the weights and 
biases. For each output unit (Y k , k = 1 m), updates its bias and weights (/ = 1 p): 

w *= %,.(° ld ) + Aw >* 

For each hidden unit (Z y , j = 1 p), updates its bias and weights (/ = 1 n): 

v<j= ^-(old) + A Vq 

One iteration through all the records in a dataset is called an epoch. Squared error e = (f* - y k ) (t k - 
y k ) can be calculated for each epoch or after a few epochs. The stop condition can be set based 
on the maximum epochs, an error threshold, or when the squared error starts increasing in either 
the training dataset or testing dataset if applicable. Once the net is trained and biases and 
weights are obtained, the feedforward algorithm is used for prediction. 

3.3 Activation function 

According to Fausett (1994), an activation function for a backpropagation net should contain 
several important characteristics: It should be continuous, differentiable, and monotonically 
nondecreasing. Furthermore, for computational efficiency, it is desirable that its derivative be 
easy to compute. For the most commonly used activation functions, the value of the derivative 
can be expressed in terms of the value of the function. Usually, the function is expected to 
saturate, i.e., approach finite maximum and minimum values asymptotically. 

One of the most typical activation functions is tire binary sigmoid function, which has range of 
(0,1) and is defined as 


/(*) = 


with 


/'(x)=/(x)[1-/(x)]. 

The binary sigmoid can have its range expanded and shifted so that it maps the real numbers into 
the interval [a, b] for any a and b. To do so, given an interval [a, b], we define the parameters 

y = b - a, 

H = -a. 


Then the sigmoid function 


9 (x) = Y / (x) - n 

has the desired property, namely, its range is (a, b). Furthermore, its derivative also can be 
expressed conveniently in terms of the function value as 

g' (x) =_L [q + g (x)][y - n - 9 (x)] 

Y 

The steepness of the logistic sigmoid can be modified by a slope parameter o. This more general 
sigmoid function (with range between 0 and 1) is 

/(x)= ! — . 

1 + exp(-for) 
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with 


/'(x) = o/(x)[1-/(x)]. 


or 


sM = Y/W-n 


and 


1 + exp(-ox) 

g'(x)=_£ L [n + g(x)iY-n -g(4 


4 Operational model design 

A neural net was designed to simulate and predict the probabilities of land transformation from 
the rural state to urban use. Land use was classified into only two categories: urban and non- 
urban, coded 1 and 0 respectively. Therefore, only one output unit is needed. Figure 2 depicts the 
backpropagation neural network land use model tailored to the unique environment of the Myrtle 
Beach region of South Carolina. 


Reliel/Slope* 
Water/Wetland* 
Distance to Major Roads 
Distance to Local Roads 
Distance to Nodes 
Distance to Water Lines 
Distance to Sewer Lines 
Distance to Urban Area 
Cost Distance to CBD 
Corporate Area 
Population Density 
Median House Value 
Distance to Waterfront 
Protected Lands* 



Neural Network 
Land Use Model 

Myrtle Beach Region 


Notes: a. Not a factor in the flat Atlantic 
coastal region. 

b. Not included in this model but 
used in the rule-based model. 


Land Use 


• Network Design: 1 lxl lxl 

• Model: Standard Backpropagation 
•Activation Function: Sigmoid 

/ (x) = X / [1 + exp ( x) ] 

/'(*) = /CO [!-/(*)] 

Range (0,1) 

•GUI: VB A/ ArcGIS/Arc View/Excel 


Inp ut Units Hi dden Units 


Output Unit 


Figure 2. Neural network land use model for the Myrtle Beach region of South Carolina, USA. 
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The input units are the site attributes and environmental factors that affect land use decision 
making and spatial distribution (Batty and Xie, 1994; Wu and Webster, 1998, Landis and Zhang, 
1997; Allen & Lu, 2003). They are grouped into five categories: market factors, physical 
suitability, service accessibility, initial conditions, and policy constraints. They represent the key 
aspects of social, economic, and environmental factors that affect land use and urban growth. In 
this study, a total of 14 predictor variables were used for building the logistic regression model 
and the neural network. Note that relief or slope is not a significant factor that causes the 
differentiation between land uses in this flat area on the Atlantic coast. It is also logical to exclude 
water features, wetlands, and protected lands first before generating sample datasets. Therefore, 
only 1 1 variables were actually used for analysis and the neural net has 1 1 active input units. 

Units of the hidden layer can be set during the training process. For simplicity, 1 1 hidden units 
were used in the final model. To make the results comparable with those of the logistic regression 
model, binary sigmoid function was used as the activation function so that the output from the 
neural net would fall between 0 and 1 . This ran be treated as a land transition probability that is 
useful for simulation and prediction of future urban growth. The cutoff value for classification is 
0.5, which is similar to that used in the logistic regression model. The cutoff values for simulating 
future urban growth were calculated based on the growth ratio as used in the previous two 
studies (Allen & Lu, 2003). 

The graphic user interface of the neural network model was developed using Visual Basic 6.00 
and runs on a PC. It was integrated with a land use module developed in ArcGIS/ArcView. GIS 
was used to generate training and testing samples and export them to the neural network model. 
The bias and weights obtained through the training or learning process were imported back to the 
GIS for a grid-based feedforward prediction so that the result can be mapped out or visualized. 
The operational neural network module provides the flexibility for users to select certain variables 
and set the corresponding units as desired. The output from the module is also in the Excel 
format and can be readily put in a word document, facilitating report and paper writing. Samples 
of the Graphic User Interface of the module are shown in Appendix A. 

5 Stratified random sampling 

Two sets of spatial data were prepared for the two baseline years, 1 990 and 2000. The 2000 land 
use layer, which has only two classes, urban (1) and nonurban (0), was used as the target layer 
in the neural network model and dependent variable in the logistic model. Spatial data for 1990 
were used for deriving datasets for the independent variables of the logistic model and the input 
units of the neural network model. The source data are shown in Table 1. Each variable grid was 
scaled to a range between 0-1 using the technique developed by Gong (1996). However, the 
range and minimum of each variable grid were stored in a separate file so that they could be 
retrieved and used for scaling the new set of grids for predicting the future growth. All of these 
data preparation processes were conducted within the ArcGIS (ESRI). 

It is neither necessary nor desirable to use the entire dataset for training the neural networks. The 
number of training patterns needed, P, depends on the number of the weights to be trained, W, 
and the accuracy of classification expected, e. According to Fausett (1994), sample size or 
enough training patterns P can be determined by the condition 

W/P = e. 
or 

P=W/e. 


For a net with w = 144 and e = 0.1 , 1440 cases or training patterns are required to be assured of 
classifying 90% of the testing patterns correctly, assuming that the net was trained to classify 
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95% of the training patterns correctly. With e = 0.05, the required training patterns are 2880. Due 
to spatial differentiation in urban growth, a sample size of 4000 cases was selected for both 
training and testing datasets. 

A stratified random sampling method was used to extract subsets from the complete coverage. 
This method guarantees that the resultant samples represent different land use classes, different 
urban patterns, or different levels of urbanization in the region. There are many ways to 
accomplish this task. In this study, a large base sample set was created then was split into nine 
subsets in two different ways for two different purposes. 


Table 1 . Operational predictor variables for the logistic model and network model. 


Variable Name 

Major Data Source 


Market factors 
Population density 

Census Blockgroup (1990, 2000) 


Housing unit value 

Census Blockgroup (1990, 2000) 


Physical suitability 
Water and wetland 

NWI (1989) and Land Cover (1990 and 2000) 


Slope /relief * 

DEM (USGS) 


Distance to waterfront 

NWI (1989) 


Service accessibility 
Distance to primary roads 

Roads (USGS 1990, 2000) 


Distance to local roads 

Roads (USGS 1990, 2000) 


Distance to major nodes 

Roads (USGS 1990, 2000) 


Distance to waterline 

South Carolina Chambers of Commerce (1990) 


Distance to sewage 
Cost distance to CBD 

Multiple sources including roads (USGS 1990) 


Initial conditions 


- 

Distance to existing urban 

Land Cover (1990, 2000) 


Policy constraints 
Protected lands 

Gap Analysis (USGS, 2003) (GA, NC, and SC) 


Corporate area 

Tiger/Lines (1990) 



First, a cell-based technique was developed to randomly select individual cells from the 2000 land 
use grid of the region. With the help of an Avenue Script, the Grid.MakeRandom function of the 
Spatial Analyst of ArcGIS was used to generate a grid with uniformly distributed random values 
between 0 and 1 . A subset of this grid, masked by the area of interest, was iteratively extracted 
either in the descending order or ascending order until a desired sample size (in number of cells 
or in percentage 25%) was obtained. This grid was reclassified to create a random mask grid. 
Then, the Grid.Sample function was used to convert all the variable grids (Table 1 ) into an ASCII 
format that can be used directly for the logistical analysis and network training. 
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From the above base sample dataset, which has 36,500 cases, nine subsets were created using 
a stratified method. Records of the dataset were numbered repeatedly from 1 to 9 and the 
records with the same number were extracted to form a subset with approximately 4072 cases. 
Each subset is capable of representing the whole population or urban reality in the region, but 
multiple datasets make it possible to use some for training and some for testing and to examine 
the global effect or strength of regional models. These subsets are called global samples. 

The second nine subsets are called areal samples or local samples. They were extracted from 
the base set following the record order and each has 4000 cases. Because the records start from 
the upper left comer to the lower right comer of the region when a grid is converted to an ASCII 
file, the resultant datasets in essence are area-based, representing different spatial patterns of 
urban space in the region. They are, therefore, called areal samples. Based on the level of 
urbanization, areas that these subsets represent can be roughly grouped into four categories: (a) 
remote rural areas with scattered urban cells (<10%, 3 subsets); (b) less urbanized areas with 
small towns (10-20%, 2 subsets); (c) urban fringe areas around metropolitan centers (20-30%, 1 
subset); and (d) urbanized areas with large green space (wetlands) along the coast (> 30%, 2 
subsets). Although these classifications are quite arbitrary, they are useful for testing the reliability 
and validity of both models with a perspective of spatial variation. 

6 Model training process 

There are a total of 18 datasets in two categories, each category having 9 sets and each set 
having approximately 4000 cases. Both logistic regression and neural network were calibrated 
and trained with the each dataset while other sets in the same category were used for testing the 
validity of each model for different geographical areas, different spatial patterns, and different 
samples. 

It was found that for each dataset, the logistic model was capable of generating the same 
regression coefficients and the same urban classifications regardless the number of runs or 
iterations. However, the results (network weights and classification accuracy) of the neural 
network model varied substantially with different runs and different parameters (such as error 
threshold, maximum epochs, and learning rates) for the same dataset. To be consistent, all 
trainings of the neural network started with maximum epochs of 3000, an error threshold of 0.05, 
and a learning rate of 0.2. 

Five strategies were used to determine which results were to be used. First, if the prediction 
results were equal or better than those from the logistic model when the maximum epochs were 
reached and the training was stopped, the weights and classifications were used as the indicator 
of the power of the model. This was the case for 1 4 of 1 8 datasets. Second, when the error 
threshold was reached quickly and the training stopped as in 4 cases of global datasets, the 
results were saved but the model was retrained with an error threshold of 0.01 . No training has 
reached that level of accuracy. Third, when the neural network generated a lower accuracy than_ 
the logistic model did for the first run, either maximum epochs were increased or the error 
threshold was lowered for continuous training to explore the potential of the neural network. No 
more than three additional runs were needed. Fourth, when the training failed in the first run as 
happened to two area-based datasets (area-sample 1 and area-sample 3) due to improper initial 
random weight or too few urban cases in the samples, additional runs were carried out. For area- 
sample 1 , desired results were obtained after resetting the initials (or neural network parameters). 
For area-sample 3, the training became successful after the learning rate was reduced to 0.05. 
Finally, when the training failed in the first few runs and the technique mentioned above did not 
work, the dataset was split into two smaller sets of odd and even numbers. Each subset was 
used to train the neural network and the resultant weights were used as the initial weights for 
training the model again using the whole dataset. This was the case for area-based dataset 5. 
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7 Prediction success rates 

7.1 Global effects 


Both models were calibrated and tested using 9 sets of sample data representing the whole 
population of the region to examine to their global effects in terms of prediction success rates. 
Results are summarized in Table 2. As these sample datasets indicate, urban area in this region 
covers about 21 .86-23.33% of the total land area. Both models generated fairly good overall 
predictions and nonurban predictions, with 88.21-90.01 % and 94.46-95.88% respectively for the 
logistic model and 89.75-93.59% and 90.13-98.37% respectively for the neural network model 
(Figure 3 and Figure 4). However, the neural networks proved to be superior to the logistic 
regression model in all but one case. The largest improvement takes place in the urban category 
with an average 4.73 net percentage. This number becomes notably larger at 10.18 if sample 
dataset 7 is excluded. More importantly, the neural network has exhibited its capability to 
differentiate urban land use in the dominantly rural environments (Figures 3 through 6). 


Table 2. Comparison of prediction success rates between the neural networks and the 
logistic regression model as trained using the global sample datasets. 


Data 


Set Land Use N (Vo) Urban 


Predicted 


Neural Networks 

Non- Correct 

urban (%) 


Logistic Regression 

Non- Correct 
urban (%) 


Urban 


Net 

Change 


1 

Urban 

939 23.05 

685 

254 

72.95 

631 

308 

67.20 

5.75 


Non-urban 

3134 76.95 

73 

3061 

97.67 

148 

2986 

95.28 

2.39 



4073 



91.97 



88.80 

3.17 

2 

Urban 

917 22.51 

760 

157 

82.88 

649 

268 

70.77 

12.11 


Non-urban 

3156 77.49 

104 

3052 

96.70 

139 

3017 

95.60 

1.10 



4073 



93.59 



90.01 

3.58 

3 a 

Urban 

925 22.72 

699 

226 

75.57 

640 

285 

69.19 

6.38 


Nonurban 

3147 77.28 

87 

3060 

97.24 

144 

3003 

95.42 

1.82 



4072 



92.31 



89.46 

2.85 

4 

Urban 

890 21.86 

675 

215 

75.84 

577 

313 

64.83 

11.01 


Nonurban 

3182 78.14 

61 

3121 

98.08 

131 

3051 

95.88 

2.20 



4072 



93.22 



89.10 

4.12 

5 

Urban 

940 23.08 

710 

230 

75.53 

646 

294 

68.72 

6.81 


Nonurban 

3132 76.92 

105 

3027 

96.65 

162 

2970 

94.83 

1.82 



4072 



91.77 



88.80 

2.97 

6 a 

Urban 

949 23.31 

724 

225 

76.29 

656 

293 

69.13 

7.16. 


Nonurban 

3123 76.69 

63 

3060 

97.98 

173 

2950 

94.46 

3.52 



4072 



92.93 



88.56 

4.37 

7 

Urban 

950 23.33 

647 

303 

68.11 

646 

304 

68.00 

0.11 


Nonurban 

3122 76.67 

51 

3071 

98.37 

153 

2969 

95.10 

3.27 



4072 



91.31 



88.78 

2.53 

8 

Urban 

928 22.79 

1326 

162 

89.11 

614 

314 

66.16 

22.95 


Nonurban 

3144 77.21 

248 

2264 

90.13 

166 

2978 

94.72 

-4.59 



4072 



89.75 



88.21 

1.54 

9 

Urban 

932 22.89 

702 

230 

75.32 

616 

316 

66.09 

9.23 


Nonurban 

3140 77.11 

97 

3043 

96.91 

143 

2997 

95.45 

1.46 



4072 



91.97 



88.73 

3.24 
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Note: (a) Training was successful after several trials; (b) Sample set 5 was first split into two 
subsets using one set for initial training. The resultant weights were used as initial training for the 
whole set. 



Figure 3. Overall success rates of both urban and nonurban classifications as predicted 
by the neural networks and logistic regression using the global datasets. 



Figure 4. Accuracy of classifications of nonurban use as predicted by the neural networks 
and logistic regression using the global datasets. 
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Figure 5. Difference of classification accuracy between the neural networks and logistic 
regression against global datasets. The trend shows that as the area becomes more 
urbanized, the difference between the two models becomes smaller. 



Figure 6. Difference of classification accuracy between the neural networks and logistic 
regression against global datasets. The trend shows that as the area becomes more 
urbanized, the difference between the two models becomes smaller. 
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The logistic model generated a slightly better prediction for the non-urban use (94.72% or 4.59 
net change) at the cost of urban land use (66.16% or -22.95 net change). Although the neural 
networks were retrained with better predictions (urban 75.1 1 %, nonurban 96.22%, and overall 
91 .40%) than those generated from the logistic model, the results of the first run were kept in the 
table because it raises an issue of neural network training to be discussed later. 

Although prediction success rates vary substantially across different land use categories, the 
accuracy of prediction is similar across different datasets. This is also true with the testing results. 
Table 5 shows the classification accuracy for each dataset based on the model calibrated using 
the first dataset. The maximum standard deviation in prediction success rates (1 .441 ) was found 
in the urban category with the smallest in the nonurban category. This was expected because the 
nonurban use is still the most dominant land use patterns in the region. It was found that for each 
testing global dataset, the prediction success rates in each category are slightly smaller than that 
for the training dataset, suggesting the neural network model best fit the dataset used for its 
training. Whether this is also true to the logistic model remained untested. 


Table 3. Prediction success rates of the neural network models trained by sample dataset 
1 and tested by eight other sample datasets. 


Dataset 

Urban 

Cells 

Nonurban 

Cells 

Total 

Cases 

Urban 

Correctly Classified (%) 
Nonurban 

Overall 

1 

939 

3134 

4073 

72.95 

97.67 

91.97 

2 

917 

3156 

4073 

71.76 

96.83 

91.19 

3 

925 

3147 

4072 

70.05 

96.28 

90.32 

4 

890 

3182 

4072 

70.45 

95.85 

90.30 

5 

940 

3132 

4072 

69.15 

95.95 

89.76 

6 

949 

3123 

4072 

70.39 

96.25 

90.23 

7 

950 

3122 

4072 

69.68 

96.09 

89.93 

8 

928 

3144 

4072 

68.43 

95.71 

89.49 

9 

932 

3140 

4072 

69.10 

96.18 

89.98 

Mean 




70.35 

96.36 

90.43 

Min 




68.43 

95.71 

89.49 

Max 




72.95 

97.67 

91.97 

Standard 




1.441 

0.626 

0.788 

Deviation 






- 


7.2 Spatial differentiations 

Prediction success rates or classification accuracy was used as the measure of power or 
reliability of the two models of concern. For each area-based dataset, the overall success rate 
and the two categorical success rates were calculated for each model and the results were 
summarized in Table 4. In general, both models have obtained good or very good overall 
predictions for the nonurban land use with an average accuracy between 89.41 % and 94.74%. 
Since the rural land covers more than a quarter of the total area in the region, prediction accuracy 
for non urban land use in terms of percentage correctly classified is often higher than for any 
other uses. Similarly, the difference in effectiveness between these two models is only marginal in 
both nonurban and all categories even though the average success rates of the neural network 
model are better than those of the logistic model (by 0.64 and 2.35 percentage points 
respectively. 
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Table 4. Prediction success rates: binary logistic regression vs. backpropagation neural 

network using area-based samples. 


Data 

Set Land Use N 


Predicted 


Urban 


Neural Network 
Non- Correct 
urban (%) 


Logistic Regression 
Non- 
urban 


Urban 


Net 

Correc Chan 9 e 
t (%) 


1 a 

Urban 

152 3.73 

92 

60 

60.53 

43 

109 

28.29 

32.24 


Non-urban 

3919 96.27 

5 

3914 

99.87 

9 

3910 

99.77 

0.10 



4071 



98.40 



97.10 

1.30 

2 

Urban 

569 14.23 

472 

97 

82.95 

396 

173 

69.60 

13.35 


Non-urban 

3431 85.78 

81 

3350 

97.64 

52 

3379 

98.48 

-0.84 



4000 



95.55 



94.38 

1.18 

3 

Urban 

267 6.68 

158 

109 

59.18 

107 

160 

40.07 

19.11 


Nonurban 

3733 93.33 

13 

3720 

99.65 

2 

3731 

99.95 

-0.30 



4000 


97.15 

96.95 



95.95 

1.00 

4 

Urban 

470 11.75 

329 

141 

70.00 

80 

390 

17.02 

52.98 


Nonurban 

3530 88.25 

48 

3482 

98.64 

27 

3503 

99.24 

-0.60 



4000 



95.28 



89.58 

5.71 

5 b 

Urban 

324 8.10 

164 

160 

50.62 

66 

258 

20.37 

30.25 


Nonurban 

3676 91.90 

26 

3650 

99.29 

16 

3660 

99.56 

-0.27 



4000 



95.35 



93.15 

2.20 

6 

Urban 

1181 29.53 

770 

411 

65.20 

786 

395 

66.55 

-1.35 


Nonurban 

2819 70.48 

51 

2768 

98.19 

258 

2561 

90.85 

7.34 



4000 



88.45 



83.68 

4.78 

7 

Urban 

1587 39.68 

1401 

186 

88.28 

1290 

297 

81.29 

6.99 


Nonurban 

2413 60.33 

551 

1862 

77.17 

350 

2063 

85.50 

-8.33 



4000 



81.58 



83.83 

-2.25 

8 

Urban 

1488 37.20 

1326 

162 

89.11 

1253 

235 

84.21 

4.90 


Nonurban 

2512 62.80 

248 

2264 

90.13 

339 

2173 

86.50 

3.63 



4000 



89.75 



85.65 

4.10 

9 

Urban 

2337 50.18 

2008 

329 

85.92 

1980 

357 

84.72 

1.20 


Nonurban 

2320 49.82 

392 

1928 

83.10 

508 

1812 

78.10 

5.00 



4657 



84.52 



81.43 

3.09 


Note: (a) Training was successful after several trials; (b) Sample set 5 was first split into two 
subsets and using one set for initial training. The resultant weights were used as initial training for 
the whole set. 


As shown in Figures 7 through 9, the most remarkable improvement was found in the urban use 
category. For the nine area-based datasets, the average prediction accuracy is 74.42% for the 
neural network model, compared to 54.68% for the logistic model, which is an improvement of 
17.74 percentage points. The largest improvement is by 52.98 percentage points in set number 
four. This suggests that the neural network model is superior in its ability to differentiate the urban 
land use from the nonurban use. On the other hand, the geographic variation of prediction 
success rates is relatively lower for the neural net (50.62-89.1 1 %) than for the logistic model 
(17.07-84.72%). No accuracy is found 
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Figure 7. Overall success rates for both urban and nonurban classifications as predictedby 
the neural networks and logistic regression using area-based datasets, (a) The datasets 
were ranked based on the level of urbanization and thus the series numbers are not the 
same as the sample numbers; (b) overall prediction accuracy tends to decline as the area 
is more urbanized; and (c) the improvement in prediction is only marginal. 


below 50% despite of great variation in the level of urbanization. All these indicate that the neural 
network-based prediction is also more stable and thus more reliable. In addition, the findings that 
success rates vary with different areas suggest that sub models or multiple models are needed 
for better approximation of a regional urban system. Existing land use models do not have the 
capability to apply different rules for addressing spatially differentiated patterns and temporally 
varying growth rates. 

It is worth reporting that the classification accuracy also changes with training sample size. During 
the preliminary study, the base dataset was split into seven subsets with the largest one having 
6657 cases while each of the rest have 5000 cases. This is the only dataset with which the 
neural network model was less successful in predicting urban use (-6.00) as shown in Table 5. Is 
this always true in more urbanized areas, or is it because the network was not adequately trained, 
or something else? To answer these questions, several other runs were performed. Some better 
results were obtained, but the accuracy was always smaller than that generated by the logistic 
model. The dataset was then broken into two smaller sub sets for retraining the net. To our 
surprise, the neural network model generated much better predictions in every category for each 
smaller subset than the logistic model did (Table 5). This appears to suggest that the smaller the 
sample size, the better the prediction based on the neural network. To further confirm this finding, 
other original datasets were subdivided and trained. Similar trends hold even though the 
difference varies from one case to another. 
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Figure 8. Success rates for nonurban use classifications as predicted by the neural networks and 
logistic regression using area-based datasets, (a) Both models have difficulties in differentiating 
the nonurban use form the urban use; (b) There is not significant difference in prediction 

accuracy between the two models. 



Figure 9. Success rates of urban use classifications as predicted by the neural networks 
and logistic regression using area-based datasets, (a) The datasets were ranked based on 
the level of urbanization and thus the series numbers are not the same as the sample 
numbers; (b) the improvement of prediction becomes greater as the area represented by 

the sample become more rural. 


Table 5 Prediction Success Rates: Binary Logistic Regression vs. Backpropagation Neural 

Network, Based on Split Smaller Samples 


Data 


Set Land Use N (%) Urban 


Predicted 


Logistic Regression 
Non- Correct 
urban 


Urban 


Neural Network 
Non- 
urban 


Correct (%) 


Net 
Change 


7 

Urban 

3176 

47.71 

2751 

425 

86.60 

2560 

616 

80.60 

(-6.00) 


Nonurban 

3481 

52.29 

681 

2800 

80.40 

369 

3112 

89.40 

9.00 



6657 




83.40 



85.20 

1.80 

7a 

Urban 

1658 

55.29 

1165 

176 

86.88 

1232 

109 

91.87 

4.99 


Non-urban 

1341 

44.71 

239 

1419 

85.59 

152 

1506 

90.83 

5.24 



2999 




86.16 



91.30 

5.14 

7b 

Urban 

1823 

49.84 

1565 

270 

85.29 

1597 

238 

87.03 

1.74 


Non-urban 

1835 

50.16 

428 

1395 

76.52 

267 

1556 

85.35 

8.83 



3658 




80.92 



86.19 

5.27 


7.3 Predicted future growth 

The projected future population growth and urban area growth in Georgetown and Horry Counties 
is summarized in Table 6 through Table 9. Although the growth ratios in both counties are not as 
large as those in the Charleston region, they represent two unique growth patterns in coastal 
South Carolina. Georgetown County, which has strict zoning ordinances, has experienced little 
growth in both population and area over the last decade (Table 6 and Table 7). In fact, 
Georgetown's population increase over the past century, only ranks sixth among eight coastal 
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counties in South Carolina. As a result, it has become the second smallest county in population in 
the coastal region. It is projected that its population will increase by only 3000 in the next thirty 
years. The urban area will expand by less than 5 square miles over the same period. On the 
contrary, Horry County exhibits an uncontrolled, sprawl growth pattern. It is the fastest growing 
county in the State of South Carolina. In 1900, Horry County’s population was almost the same 
as Georgetown’s population. It has increase 8.42 times by 2000. This growth rate ranks the first 
in the entire state. Over the last decade, its urban area increased by nearly 40 square kilometers. 
By year 2030, its urban area will expand from 142.84 km 2 in 2000 to 259.21 km 2 at the current 
pace and more likely triple to 484.65 km 2 at the 3:1 growth ratio (Table 8). 

Figure 9 and Figure 10 show the spatial extent of historical urban growth during 1990-2000 and 
predicted urban growth through 2030. In comparison, predicted urban growth for the other two 
coastal regions is shown in Figure 12 and 13 in the Appendix. Urban sprawl in this region mainly 
took a zonal pattern expanding both ways along the coast and spreading toward inland towns 
along the major highways. This spatial trend is predicted to continue through 2030. Areas around 
the major road intersections along the US 17 Bypass will continue to be the hot growth spots in 
the region. Due to the constraint of wetland areas, no significant growth is anticipated in the 
stretch between the Myrtle Beach and North Myrtle Beach, but the upland or dry land areas along 
two major roads (US Hwy 501 and SC Hwy 9) leading inland will see significant urban 
developments. The triangle area, formed by US Hwy 501 , US Hwy 17, SC HWY 707 and SC 
HWY 544, and SC Hwy 9, is the only large tract of physically developable land adjacent to the 
Myrtle Beach and thus is the prime target for future development. This will create significant 
pressure on the coastal ecosystems in that region. Summary tables about the future urban growth 
are shown in tables six through nine. 

Table 6. Projected urban area growth in square kilometers, Georgetown County, South 

Carolina, 1990-2030. 


Year 


Population 

Growth 


Urban Area Growth in Square Kilometers 
Ratio 0.48:1 Ratio 1:1 Ratio 2:1 Ratio 3:1 


1990 

46302 

22.11 

22.11 

22.11 

22.11 

1995 

51050 

23.03 

23.03 

23.03 

23.03 

2000 

55797 

23.95 

25.17 

27.32 

29.46 

2005 

60544 

24.87 

27.32 

31.60 

35.88 

2010 

65291 

25.79 

29.46 

35.88 

42.31 

2015 

70038 

26.70 

31.60 

40.17 

48.73' 

2020 

74785 

27.62 

33.74 

44.45 

55.16 

2025 

79532 

28.54 

35.88 

48.73 

61.58 

2030 

84279 

29.46 

38.03 

53.02 

68.01 


Note: (a) Urban areas are calculated based on the projected population growth and the growth 
ratios; (b) urban area in the county will double in the next thirty years based on the 2:1 growth 
scenario. 


Table 7. Projected urban area growth in percentage, Georgetown County, South Carolina, 

1990-2030. 


Year 


Population 
Growth (%) 


Urban Area Growth in Percentage (%) 

Ratio 0.48:1 Ratio 1:1 Ratio 2:1 Ratio 3:1 
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1990 


1995 

10.25 

4.17 

4.17 

4.17 

4.17 

2000 

20.51 

8.33 

13.86 

23.55 

33.23 

2005 

30.76 

12.48 

23.55 

42.92 

62.29 

2010 

41.01 

16.63 

33.23 

62.29 

91.36 

2015 

51.26 

20.78 

42.92 

81.67 

120.42 

2020 

61.52 

24.93 

52.61 

101.04 

149.48 

2025 

71.77 

29.09 

62.29 

120.42 

178.54 

2030 

82.02 

33.24 

71.98 

139.79 

207.60 


Note: (a) Year 1990 was used as the starting baseline year for percentage calculation; (b) the 
ratio of urban area growth (%) to population growth (%) was 0.48: 1 for the period between 1990 
and 2000. 


Table 8. Projected urban area growth in square kilometers, Horry County, South Carolina, 

1990-2030. 


Population 

Growth 


Urban Area Growth in Square Kilometers 
Ratio 1 :1 Ratio 2:1 Ratio 3:1 


1990 

144053 

104.06 

104.06 

104.06 

1995 

170341 

123.45 

123.45 

123.45 

2000 

196629 

142.84 

142.84 

142.84 

2005 

222917 

162.24 

180.82 

199.81 

2010 

249205 

181.63 

218.80 

256.78 

2015 

275493 

201.03 

256.78 

313.75 

2020 

301781 

220.42 

294.76 

370.71 

2025 

328069 

239.82 

332.73 

427.68 

2030 

354357 

259.21 

370.71 

484.65 


Note: (a) Urban areas are calculated based on the projected population and growth ratios; (b) 
urban area in the county will triple in the next thirty years based on the 3:1 growth scenario. 
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Table 9. Projected urban area growth in percentage, Horry County, South Carolina, 1990- 

2030. 


Year 


Population 

Growth 


Urban Area Growth in Square Kilometers 
Ratio 1:1 Ratio 2:1 Ratio 3:1 


1990 

1995 

18.25 

18.63 

18.63 

18.63 

2000 

36.50 

37.27 

37.27 

37.27 

2005 

54.75 

55.91 

73.77 

92.01 

2010 

73.00 

74.54 

110.26 

146.76 

2015 

91.24 

93.19 

146.76 

201.51 

2020 

109.49 

111.82 

183.26 

256.25 

2025 

127.74 

130.46 

219.75 

310.99 

2030 

145.99 

149.10 

256.25 

365.74 


Note: (a) Urban areas are calculated based on the projected population growth and the growth 
ratios; (b) urban area in the county will double in the next thirty years based on the 2:1 growth 
scenario. 
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Figure 10. Urban growth in the Georgetown-Horry region of South Carolina, 1990-2000 
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Figure 11. Predicted urban growth in the Georgetown-Horry region of South Carolina, 2010 
2030. Demand of future urban sizes in the Horry County was based on a 3:1 urban area- 

population growth ratio. 


7 Conclusions 

The use of an appropriate relationship model is critical for reliable prediction of future urban 
growth, identification of proper variables and mathematic functions and determination of the 
weights or coefficients are the key tasks for building such a model. Although the conventional 
logistic regression model is appropriate for handing land use problems, it appears insufficient to 
address the issue of interdependency of the predictor variables. This study used an alternative 
approach to simulation and modeling urban growth using artificial neural networks. It developed 
an operational neural network model trained using a robust backpropagation method. The model 
was applied in the Myrtle Beach region of South Carolina, and tested with both global datasets 
and areal datasets to examine the strength of both regional models and areal models. The results 
indicate that the neural network model not only has many theoretic advantages over other 
conventional mathematic models in representing the complex urban systems, but also is 
practically superior to the logistic model in its capability to predict urban growth with better 
accuracy and less variation. The neural network model is particularly effective in terms of 
successfully identifying urban patterns in the rural areas where the logistic model often falls short. 
It was also found from the area-based tests that there are significant intra-regional differentiations 
in urban growth with different rules and rates. This suggests that the global modeling approach, 
or one model for the entire region, may not be adequate for simulation of a urban growth at the 
regional scale. Future research should develop methods for identification and subdivision of these 
areas and use a set of area-based models to address the issues of multi-centered, intra- 
regionally differentiated urban growth. 
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Figure 12. Sample graphic user interfaces of the neural network land use model. 
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Appendix A. Predicted urban growth in the BCD and BCJ regions 



Figure 13. Predicted urban growth in the BCD (Berkeley, Charleston and Dorchester) and BCJ 
(Beaufort, Colleton and Jasper) regions of South Carolina, 2005-2030. 
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Figure 14. Predicted urban growth in 2030 in the BCD and BCJ region. 
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